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ABSTRACT 

In  this  paper,  which  is  part  II  in  a  series  of  two,  the  investigation  of  the  Galerkin 
finite  element  solution  to  the  Helmholtz  equation  is  continued.  While  part  I  con¬ 
tained  results  on  the  h- version  with  piecewise  linear  approximation,  the  present 
part  deals  with  approximation  spaces  of  order  p  >  1.  The  method  is  assumed  to 
be  uniform  both  w.r.  to  h  and  p.  As  in  part  I,  the  results  are  presented  on  a 
one-dimensional  model  problem  with  Dirichlet/  Robin  boundary  conditions.  In 
particular  there  are  proven  stability  estimates,  both  w.r.  to  data  of  higher  regu¬ 
larity  and  data  that  is  bounded  in  lower  norms.  The  estimates  are  shown  both 
for  the  continuous  and  the  discrete  spaces  under  consideration.  Further  there 
is  proven  a  result  on  the  phase  difference  between  the  exact  and  the  Galerkin 
finite  element  solutions  for  arbitrary  p  that  had  been  previously  conjected  from 
numerical  experiments.  These  results  and  further  preparatory  statements  are 
then  employed  to  show  error  estimates  for  the  Galerkin  finite  element  method. 
It  becomes  evident  that  the  error  estimate  for  higher  approximation  can  -  with 
certain  assumptions  on  the  data  -  be  written  in  the  same  form  as  for  the  piece- 
wise  linear  case,  namley,  as  the  sum  of  the  error  of  best  approximation  plus  a 
pollution  term  that  is  of  the  order  of  the  phase  difference.  The  paper  is  concluded 
by  a  numerical  evaluation. 
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1  Introduction 

•  xu  A  r.prf  rof  two'l  of  an  investigtion  devoted  to  the  numerical 

This  paper  is  the  second  part  (  ^  reduced  wave  (Helmholtz)  equa- 

analysis  of  the  Galerkin  finite  elemen  m  years-  related  results  have 

following  one-dimensional  model  problem. 

Let  H  =  (0, 1)  and  consider  the  boundary  value  problem, 


u’\x)  +  kM^)  = 

u(0)  =  0 

u'(l)-ifcu{l)  =  0. 


(1.1) 

(1.2) 

(1.3) 


(1.4) 


or.  equivalently,  the  variational  problem:  Find  u  €  u(0)  =  0  such  that 

B(u.v)  =  //  (u'(x)u'(x)  -  =  r  /(*)'(*)* 

''°'tintti’’w1lS-‘&«^  ‘he  cane  of  large  wavenumber  we  aaaume 

throughout  the  paper  that  fc  >  1. 

We  showed  [IB]: 

Al.  For  data  /  €  I‘(0. 1),  the  BVP  (1.1-1.3)  haa  unique  solution  u  €  K‘(0, 1). 

A2.  For  data  /  S  H‘(0, 1)',  the  VP  (1.4)  has  unique  solution  u  €  JP'(0. 1). 

A3.  The  form  B  defined  by  (1.4)  is  continuous  for  each  fixed  wavenumber  h 

\B{u,v)\  <  C'o(fc)|u|ilu|i 

with  Co  =  1  + 

A4.  The  Babuska-Bresui  constant  q  of  the  VP  (1.4)  is  inversely  proportional  to  the 

wavenumber  k. 

Remark  1=  Obviously,  the  statements  A  l^d  A2  hold 

the  data;  in  general,  the  operators  given  by  BVP  (LM-d)  and  vr  ^  ; 
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mappings  from  to 

For  the  Galerkin  finite  element  method  with  piecewise  linear  approximation  (p  =  1) 
we  then  proved: 

FI.  The  finite  element  solution  is  (asymptotically)  quasioptimal  w.r.  to  the  H^- 
seminorm  provided  Ph  is  sufficiently  small,  where  h  is  the  size  of  the  elements. 

F2.  The  discrete  B-B-condition  holds  with  ~  p 

F3.  In  the  preasymptotic  range,  the  error  bound  in  -seminorm  consists  of  a  term 
of  order  kh  (reflecting  the  error  of  best  approximation)  and  a  pollution  term 
of  order  k^h?  (reflecting  the  phase  difference  between  exact  and  finite  element 
solutions). 

Finally,  we  conluded  from  a  numerical  evaluation  that  the  terms  mentioned  in  F3 
do  indeed  occur  in  the  error  and,  consequently,  the  estimate  (obtained  from  F3.  by 
setting  6  =  kh) 

\eU<c^e  +  C2ke\  (1.5) 

where  C\B  is  the  minimal  error  in  the  approximation  space,  cannot  be  principally  im¬ 
proved. 

In  the  present  paper  we  extend  our  study  to  elements  of  arbitrary  (but  fixed)  poly¬ 
nomial  order  p.  Due  to  the  oscillatory  character  of  the  propagating  solutions,  the 
application  of  higher  order  (egs.,  quadratic  or  cubic)  elements  is  considered  a  natural 
choice  in  many  applied  computations  (cf.  [TPl],  [TP2]  and  references  therein).  How¬ 
ever,  except  for  an  asymptotic  error  estimate,  no  results  stemming  from  the  numerical 
analysis  of  Galerkin-type  FEM  for  Helmholtz  problems  seem  to  be  available  for  p  >  1. 
In  the  abovementioned  cisymptotic  estimate,  which  is  given  without  rigorous  proof  in 
[BayGTu],  it  is  assumed  that  k^h  <C  1  (cf.  Theorem  3.2  of  this  paper).  As  we  showed 
in  part  I,  this  assumption  on  the  meshsize  does  not  cover  the  size  of  practically  applied 
meshes  where  the  wavenumber  k  (but  not  which  is  a  significant  difference  for  large 
k)  is  commonly  normalized^  by  the  stepwidth  h  [TPl],  [HH]. 

Consequently,  as  in  part  I,  it  is  the  principal  goal  of  the  present  investigation  to 
show  error  estimates  that  hold  under  assumptions  on  hk  only.  We  call  these  esti¬ 
mates  "preasymptotic”  in  order  to  distinguish  the  results  from  the  abovementioned 
asymptotic  estimates.  In  section  3  of  this  paper,  we  prove  two  error  estimates  that 

^i.e.  h  is  choosen  such  that  the  product  hk  remains  close  to  some  apriori  fixed  value 
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hold  in  the  preaeymptotic  range  Both 

‘’‘n^To^feC  wtrilCfr^  ^  PoltaWo''  '>5’.  "t"® 

Li=-estimateofthedata,t^^^^^^ 

the  pollution  term  of  the  second  estimate  (in  h,  h  and  p)  in 

““^irturns  out  (both  theoretically  and  in  the  “"^X‘ot\rh^^^ 

in  the  best  approximation  estimate  grows  on  y  moderately  ?  _  theoretically 

the  constant  in  the  esBmate  of  ^n  orsl^ed  in 

?rrlli"exT^^^  “  °P“ 

is  sharp  w.r.  to  p. 

investigate  the  Galerkin  finite  element  of  the  finite 

tion  3.1)  we  identify  the  W'f subStion  3  3  we  prove  the  proposition 
element  solution  proce  3,4  jj  aevoted  to  stability  propositions  on 

on  the  phase  lag  (Theorem  3.2).  {qj.  general  L^-data 

Srt^jTas  follows’-  in  the  approximaBon  — 
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interpolation  properties  w.r.  to  u.  It  is  these  interpolation  properties  that  we  make 
use  of  in  the  error  analysis  where  we  show  that  the  finite  element  error  is  majorized 
by  certain  norms  of  the  difference  u  —  s.  Since  u  —  s  occurs  as  the  data  in  variational 
problems  related  to  the  error  estimation,  we  have  to  prepare  the  error  analysis  with 
stability  specific  propositions.  Thus  prepared,  we  turn  to  the  error  analysis  of  the 
Galerkin  finite  element  method  (subsection  3.5).  We  show  an  asymptotic  estimate  a.nd 
then,  in  light  of  the  previous  discussion,  turn  to  a  study  of  the  finite  element  solution 
in  the  preasymptotic  range. We  show  that,  for  oscillating  solutions,  the  pollution  term 
is  of  the  order  of  the  phase  lag. 

Finally,  in  section  4,  we  discuss  results  of  selected  numerical  experiments  that  have 
been  carried  out  to  evaluate  the  error  estimates  of  section  3. 


2  Analytic  solution  properties 

2.1  Notations  and  Preliminaries 

We  will  use  the  following  notations: 

•  Constants:  If  not  stated  otherwise,  all  constants  C,  D,  E, . . .  or  C,-,  where  i  is  a 
natural  number,  are  understood  to  be  generic,  independent  of  all  parameters  of 
the  given  estimate,  and  having,  in  general,  different  meanings  in  different  context. 
However,  some  specific  constants  that  are  used  repeatedly  with  the  same  meaning 
in  different  context  are  marked  by  literal  subscript;  so,  we  write  Ca  for  a  constant 
arising  from  the  principal  approximation  estimate  in  the  discrete  space,  and  so 
forth. 

•  i^(0, 1)  is  the  space  of  all  square-integrable  complex-valued  functions  equipped 

with  the  inner  product  ^ 

(u,u;):=  /  v{x)w{x)dx 
Jo 

and  the  norm  - 

\\w\\  ■= 

•  H^{0, 1)  denotes  the  Sobolev  space 

H‘’  =  {u|u  €  L'^  Ad'u  e  =  l...s} 

where  d'u  are  the  derivatives  of  order  i  in  the  distribution  sense, 
define  the  subspace 

1)  =  {u  €  H^{0,  l)15’u(0)  =  a’u(l)  =  0,  i  =  0 . . .  s  - 


As  usual,  we 
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We  will  also  work  with  subspaces  and  consisting  of  functions  with  Dirich- 
let  data  0  given  only  in  a;  =  0  or  x  =  1,  resp.  By  |u|3  :=  a  seminorin^is 

given  on  H^.  A  norm  of  the  space  H^{0, 1)  is  defined  as  ||u||a  :  I^L) 

On  HI,  Hl^  and  the  seminorm  |  •  is  a  norm  equivalent  to  ||  •  ||,. 

•  1)  =  1))'  denotes  the  dual  to  space  equipped  with  the  norm 


ll/ll-» 


sup 


(2.1) 


•  As  usual,  /h),  i  =  0,1,2,...,  denotes  the  i-th  derivative  of  /.  We  generalize 
this  notation  for  integration:  f^~'^  is  a  function  s.t.  “  fi^)- 

specifically,  we  define  for  /  €  i/^(D)  and  z  =  0, 1, . . . 

/(-<-')(x)  =  -  /'  /'-■>(()*■ 

Jx 


With  these  definitions,  the  dual  i/^-norm  of  a  L^-function  /  is  equal  to  the  L  -norm 
of  /(-'): 


Lemma  2.1  For  f  €  and  /  =  0, 1, 2, . . . 

ll/ll-/  =  ll/■'^ll 


(2.2) 


holds. 


Proof  :  Let  /  =  1  and  write  F  =  Then,  by  partial  integration, 


Ill’ll  = 


sup 

v€L2{Q) 


loFv 

ll^ll 


sup 

ii€n*(n) 


kfV 


holds  with  V  :=  Jq  v(t)dt.  Obviously  V  €  and  V  —  v.  On  the  other  hand, 

every  V  G  H^g{Cl)  can  be  represented  by  an  integral  of  a  L^-function,  hence 


This  proves  the  statement  for  /  =  1.  The  induction  to  higher  /  is  obvious  and  the  proof 
is  completed.  <d 
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Variational  forms  arising  from  the  Helmholtz  equation  are,  in  general,  indefinite. 
One  can  obtain,  however,  coercive  forms  if  the  wavenumber  is  properly  restricted  (we 
write  K  in  order  to  distinguish  this  problem  from  the  general  case  where  k  may  be 
large).  For  later  use,  we  consider  here  the  case  of  Dirichlet  boundary  conditions. 

Find  u  €  V  =  such  that  for  all  u  €  V 

Bk{u,  =  ~  Jo 


holds. 


Lemma  2.2  Letu  eV  =  Hl{Q,)  be  the  solution  to  the  VP  (2.3)  with  data  f.  Assume 
that  0  <  K  <  a  Ktt,  then 


IIHl  < 

Ih'll  < 

Ill'll  < 


TT^  — 
IT 

TT^  — 


TT 


TT^  — 


ll/ll 

ll/ll 


(2.4) 

(2.5) 

(2.6) 


Proof  :  All  inequalities  are  trivial  for  u  =  0,  hence  we  may,  without  loss  of 
generality,  assume  l|u||  >  0.  It  is  easy  to  see  that,  for  all  u  G 

Ill'll  >  7r||u||. 


Hence 

and 


Bk{u,u)  > 


Bk{u,u)  > 


TT 


/||2 


From  the  first  inequality  we  easily  conclude  eq  (2.4): 


v€L2{Q)  Ikll  ll“ll 


Equation  (2.6)  follows  similarly. 
Eq  (2.5)  follows  from  (2.4)  by 
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The  proof  is  completed.  <1 

Remark  2:  The  statement  of  the  lemma  holds  also  if  u  and  v  are  choosen  from  a 
Hilbert  subspace  14  C  V.  Indeed,  since  obviously  min„ev(|w|i/||«||)  <  niinuev).(|w|i/||«ll)5 
the  form  Bk  is  strongly  elliptic  on  the  subspace  and  the  same  arguments  apply. 


2.2  Stability  estimates  for  higher  p 

In  part  I  we  proved  stability  estimates  for  L^-  and  H  ^-data.  We  now  generalize 
these  properties  in  two  directions.  Namely,  we  first  consider  data  of  higher  (J;han  L  ) 
regularity  and  prove  that  for  /  >  2  the  solution  norm  |u|/4.i  is  bounded  by 

We  then  show  a  dual  result,  i.e.  we  bound  |u|i  by  ||/  '”||.  In  this  case  we  consider 
data  of  a  specific  type  -  ’bubble’  data  -  for  which  the  integrals  vanish  at  the  boundaries 
of  fl.  The  sense  of  this  assumption  will  become  clear  in  the  error  analysis. 

Theorem  2.1  (Continuous  stability):  Let  f  be  the  data  and  u  the  solution  of  the 
BVP  (1. 1-1.3).  Assume,  for  I  >  1,  f{x)  €  1).  Then  u  €  1)  and  the 

estimate 

(2-7) 

holds  for  a  ■positive  constant  €^{1)  <  Dl,  where  D  does  not  depend  on  k  and  1. 

Remark  3:  Except  for  the  dependence  on  k  (and  /),  the  statement  is  similar  to  the 
well  known  regularity  result  for  the  Laplace  equation  (cf.  [Sch,  p.52/53]). 

Proof  :  Let  us  first  consider  the  case  1  =  2.  We  have  to  prove  |u|3  <  Ck\\f\\i.  We 
start  from  the  Green’s  function  representation  of  u  (see  part  I  for  details): 


u(x)  = 

=  /  0(x, 
Jo 

s)f{s)ds 

(2.8) 

where 

I 

(  sin  kxe'^‘:  0  ■<  x  <  s 

^(^r. 

s)  =  T  ^ 

1 

1 

• 

(2.9) 

k 

[  sinifcse’^^;  s<x<l 

By  partial  integration. 

u{x)  = 

-  ["  H{x,s)f{s)ds 

Jo 

(2.10) 

Finite  Element  Solution  to  the  Helmholtz  Equation  II .. . 
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with 


H{x,s)  :=  J  G{x,s)ds  =  ^ 


ismkxe'^^^  +  l]  0<x<s 
cos  fcsc'*^;  s  <  X  <  1 


(2.11) 


For  any  fixed  s  (or  x,  resp.),  H{x,  s)  is  a  /f^-function  of  x  (or  s,  resp.).  In  the  boundary 
points,  H{x,0)  and  H{x,  1)  are  smooth  ((7°°)  functions.  We  now  estimate 


u(a;)|  <  |i?(a;,0)||/(0)H-|i/(a:,l)||/(l)|  +  sup|/f(x,s)lll/'l|. 


From  eq  (2.11)  we  have  directly 

Va:,s:  \H{x,s)\<^ 

and  with 

Vs  :  l/(s)|  <  n/2||/||, 

(for  the  proof  see,  e.g.,  [BKS])  we  get 

M<^(l+2V2)||/||,.  (2.12) 

For  an  estimate  of  the  derivatives  of  u,  differentiate  in  eq  (2.10)  to  obtain 

u\x)  =  [H^{x,  s)f{s)Y=l  -  f  H^{x,  s)f\s)ds  (2.13) 

JO 

and  from  this  and  differentiation  w.r.  to  x  in  (2.11), 

\u\i<l{l+2V2)\\f\\i.  (2.14) 

Similarily,  since  H  €  we  obtain  differentiating  (2.13) 

Ms  <  (1  +  2^)11/11,.  (2.15) 

Finally,  since  u  €  Fl^(ll),  the  differential  equation  u'"  +  =  /'  holds  (at  least  in  the 

weak  sense).  Hence 

jula  <  +  2^:^|u|i|/|i  +  \f\v 

and  with  (2.14)  we  obtain 


\u\l<  0^^^  +  20^/111  +  11/111 
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or,  equivalently, 

Ms  <£■  Ml/ll.  (2-16) 

which  proves  the  statement  for  /  =  2. 

For  higher  I  we  proceed  analogously.  First  we  integrate  in  eq  (2.8)  by  parts  (/  -  1) 
times: 


Jo 


(2.17) 


where  G(-^)(x,s)  =  /  with  G(°)(x,s)  :=  G(x,s)  and  appropriate  inte¬ 

gration  constants  for  continuity  at  x  =  s. 

For  fixed  x  resp.  s,  the  regularity  is  now  GG^\x,  s)  G  1).  At  the  boundaries 

we  have  again  (j'^"-'^(x,0),G^“'’^(x,  1)  G  (^““(0, 1). 

Therefore  differentiation  of  GG^'>{x,s)  w.r.  to  x  is  is  well  defined  at  most  1  times. 

Hence  for  j  =  1, . . . ,  I, 

|»W(x)|  <  |G«-"(*.0)|l/(0)l  +  |C‘'"‘’(^.l)|l/(l)l+ ■••  + 

+  o)|  +  |g<'-'+’>(*.  i)|  + 

lyo  I 

The  data  is  bounded  in  any  point  by 

Vx  G  [0, 1],  Vj  :  |/•'^(a:)|  <  <  V2||/||j+i. 

Generalizing  eq  (2.11),  the  integrals  of  the  Green’s  function  can  be  written  in  the  form 
G(-’")(x,s)  =  k-^^-^^'>(p{x,s)  +  k-^Pm-i{x,s) 

where  <f{x,  s)  is  an  oscillating  part  (with  sup  |v?|  <  1)  and  Pm-i  is  a  sum  of  polynomials 
of  degree  m  —  1  in  s  and  x,  resp.  Consequently,  Pm-i  is  bounded  on  Vt  and 

[  Ci{j,m)k~'^  +  C2{j,m)k^~”'~'^  if  j<m 

|G(J-”’)(x,5)|  <  { 

'  ^  if  j  >m 
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In  particular,  for  j  —  I  —  1  there  exists  a  constant  Cs^l)  —  niaxj,m(C'i(j,  rn)) 

s.t. 

ul'-'>(i)|  <  C3  (2  («:'-=  +  k‘-'  +  ...  +  k-')  +  r‘)  \/2  ll/ll, 

Similarly,  for  j  =  I  there  exists  C'4(/)  s.t. 

\u^^\x)\  <  C,  (2  (F-2  +  F-"  +  . . .  +  1)  +  1)  ||/||p-i. 

Hence  for  fc  >  1  and  /  >  2  there  are  constants  Csil)  and  Ce{l)  not  depending  on  k  s.t. 

<  CsF-=||/||,-.  (2-18) 

Ifl,  <  CsF-l/i-.  (2.19) 

hold.  By  their  definition,  C5  and  Ce  are  of  order  /. 

Then,  finally, 

l«l?„  =  jT’ (/<'-■) 

<  k'^\u\f_i  +  2A:^|u|/_i|/|/_i  + 

<  +  2  eCsk‘-^  +  1)  ||/||?_i 

and  the  statement  readily  follows.  The  proof  is  completed.  <d 

We  now  proceed  to  a  second  stability  result  that  is  dual  to  the  first  one,  i.e.  we  will 
bound  lower  solution  norms  by  negative  norms  of  the  data.  We  employ  two  auxiliary 
problems: 

1.  Consider  the  2nd-order  Dirichlet  BVP  on  fl  =  (0, 1). 

w"  -  k^w  =  -g  (2.20) 

it;(0)  =  n>(l)  =  0  (2.21) 

and  the  associated  variational  problem:  for  g  €  H^{Q,  1)'  find  u  €  ifo(0,l)  s.t. 
Vt;e/fj(0,l):  Bi{w,v)  =  w'v' +  k"^  wv  =  gv  =  {g,v)  (2.22) 
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2.  Consider  the  4th-order  Dirichlet  BVP 

w""  +  k^w  =  g  (2.23) 

iy(0)  =  iy(l)  =  iw^(O)  =  =  0  (2.24) 

and  the  associated  variational  problem;  for  g  €  H^iO,  1)'  find  u  €  (0,1)  s-t- 

Vn  €  H^{0, 1)  :  B2{w,  v)  =  w"v"  +  ^  i  ^ 


The  forms  Bi  and  B2  are  coercive,  hence  eqs  (2.22)  and  (2.25)  have  unique  solutions. 
Lemma  2.3  Let  ui  and  U2  be  the  solutions  of  Bi{w,v)  =  {g,v)  and  B2{w,v)  =  {g,v), 


spectively. 

Then  for  ui 

llx'.ll  <  lls<-‘>ll 

(2.26) 

and  for  U2 

llu;ii  <  ii5<-">ii 

(2.27) 

(2.28) 

kV2\\  <  5ll5'-'’ll 

(2.29) 

hold. 

Proof  :  Let  |||  •  |||i  and  HI  •  III2  be  the  energy  norms  induced  by  the  forms  Bi  and 
B2,  resp.  In  these  norms,  the  B-B-constant  is  trivially  7  =  1  and  hence 


ilKI||2<|lMll2- 

It  is  then  easy  to  see  from  the  definitions  of  the  various  norms  that 

IKil<||Kllli<llM|i;<  w  = 
hold.  This  proves  (2.26)  and  (2.27). 


Finite  Element  Solution  to  the  Helmholtz  Equation  II .. . 


13 


To  show  eq  (2.29),  let  u  G  From 

0  <  (h'll  -  k^uWY 


we  conclude 

2e\\u"\\\\u\\  <  B,{u,u) 

and,  by  Schwarz  inequality  and  partial  integration, 

B2{u,u)  >  2fc^l(«,u")l  =  2P||u'|p.  (2.30) 

On  the  other  hand,  also  by  partial  integration  and  Schwarz  inequality, 

B2{u,u)  <  |(5^■'^u')l  <  lk^"'^||||«'ll- 


Hence 


2k^ur<y-^wi 


Cancelling  l|u'||  we  obtain  (2.29). 

Finally,  multiplying  (2.30)  by  the  trivial  relation  B2{u,u)  >  Iju  ||  we  have 


(2.31) 


But  is  majorized  by  and  (2.28)  readily  follows  cancelling  ||n"f 

and  taking  roots.  This  completes  the  proof.  < 


Theorem  2.2  (A  dual  stability  result):  i4ssume  that,  for  integer  m  >  1,  there  is 
given  a  function  f  €  L^Q)  such  that  /(-’)(0)  =  /("‘^(l)  =  0  for  i  =  Let 

u  €  Hl^{Q.)  he  the  solution  to  the  VP  (1.4)  with  data  f. 

Then  ,  ^  ^ 

lull  <  C'iA:'"||/(-’")||  +  C'2||/-'^||  (2.32) 

holds  with  Ci,C2  independent  of  k. 

Remark  4:  The  assumption  on  the  data  means  that  m  integrals  of  /  vanish  in  a:  =  0 
(note  that  all  integrals  vanish  at  the  endpoint  by  definition).  Without  this  assumption 
we  can  just  prove  that 

|n|.  <  C  k  ||/<-‘>|| 

(see  part  I).  In  general,  |u|i  cannot  be  bounded  by  a  term  C||/l  indepenently  of  k. 
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Proof  :  (1)  For  m  =  1  the  statement  simplifies  to 

Nl.  < 

and  is  directly  obtained  from  the  B-B-condition  (Introduction,  A4.)  and  general  theory 
([BA]). 

(2)  Let  m  =  2.  The  basic  ingredience  of  the  argument  is  the  introduction  of  a  smoother 
kernel  to  the  Green’s  function  representation  u(a;)  =  (G(a:,5),/(5))  of  the  solution.  We 
define 

K{x,s)  :=  G{x,s)- H{x,s) 

where  H{x.,s)  is  the  Green’s  function  to  the  first  auxiliary  problem  (2.22).  Per  defini¬ 
tion  of  Green’s  functions,  the  equations 

Hxx{x,  s)  -  H{x,  s)  =  -6s{x)  (2.33) 

Gx.{x,s)-\-eG{x,s)  =  -S,{x)  (2.34) 

hold  (in  the  sense  of  distributions).  Hence,  after  substraction, 

Kx,  =  -k\G  +  H).  (2.35) 

Since  G,H  £  H^{n)  (as  functions  of  x  for  any  fixed  5)  it  follows  that  K  €  We 

have  constructed  an  equivalent  integral  representation  of  u  as 

u{x)  =  f  K{x,s)f{s)ds  +  f  H{x,s)f{s)  :=  ui{x)  +  U2{x).  (2.36) 

Jo  Jo 

For  estimation  of  ||u'J|  we  integrate  by  parts 

„;(x)  =  Kx{x,s)f{s)ds  =  £  Kxss{x,s)f-^\s)ds, 

(note  that  the  boundary  terms  vanish  due  to  the  specific  assumption  on  /)  hence  by 
Cauchy- Schwarz  inequality 

From  eq  (2.35)  and  the  symmetry  property  of  Green’s  functions,  /fr..  =  {G+  H)x  and 

||^„.||<t^(||Gr||  +  l|/fr||) 
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follows.  It  is  straightforward  to  show  that  HGrll  and  ||^fi|l  are  bounded  independently 
of  k  and  we  thus  have  ^ 

l»;(x)l<CPl|/<-^>||.  (2-37) 

To  estimate  U2  we  apply  eq  (2.26),  Lemma  2.2,  for  g  =  f'. 

Ih2ll  < 

which,  together  with  eq  (2.37)  yields  the  statement,  case  2. 


(3)  Consider  the  case  3  <  m  <  6.  In  analogy  to  the  previous  step,  we  introduce  a 
still  smoother  kernel  to  the  integral  representation  of  u.  Let  us  prove  the  extreme  case 
m  =  6.  We  start  from  eq  (2.36)  and  define 

L(x,s)  :=  K{x,s)  -  J(x,s) 

where  J{x,s)  is  a  Green’s  function  with  the  same  singularity  as  K.  Then 

u(x)  =  /'  L(x,  s)f{s)ds  +  /'  J(x,  s)f{s)  +  /  if (x,  s)f{s) 

Jo  ''0 


:=  Uo(a;)  +  ui(a:)  +  U2(x).  (2.38) 

To  find  the  appropriate  J,  add  eqs  (2.33)  and  (2.34): 

(G  +  HU{x,s)  +  P(G  -  H){x,s)  =  -26{x,s). 

Multiplying  this  eq  by  and  substituting  K  =  G  -  if,  eq  (2.35)  leads  to 

s)  -  k^I<{x,  s)=^2k^  6,(x,  s). 

Hence  if  J(x,  s)  is  the  Green’s  function  to  the  second  auxiliary  problem  with  data 
5f(x)  =  2  f:^/(x)  then 

Lxxxx  =  {K  -  JU.S  =  {K  +  J)- 

Since  itT,  J  €  H^{^)  it  follows  that  L  €  H\9).  Thus,  integrating  by  parts, 

l<(^)l  =  \j  L,ssssss{x,s)f^-^\s)ds 

=  k^\j{K  +  J)s..{^.s)f^-^\s)ds 

=  k^  I  J{G  +  if  ).(x,  s)f-^\s)ds  +  J  JJ^-^\s)ds'^ 
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and  therefore 

IKII  < 

where  uz  is  a  solution  to  the  second  auxiliary  problem  (2.25)  with  data  g  =  2k 
From  Lemma  2.2  and  Poincare’s  inequality, 

IK||<|u3|2<4P1|/(-®)||, 

VjpTirp 

By  definition  (2.38),  ui  is  the  solution  to  the  VP  (2.25)  with  data  g  =  2k^f  and  U2 
solves  the  VP  (2.22)  with  data  /.  Then  it  follows  from  Lemma  2.2.  that 

and 

ll-iil  <  ll/'■"ll 

so  that  finally  u  can  be  estimated  as 

|n|.  <C.t'>||/'-''l||  +  C.||/(-'>|| 

which  proves  the  statement  for  m  =  7. 

For  m  =  4  . . .  m  =  6  the  argument  is  completely  analogous,  the  only  difference  being 
the  number  of  partial  integrations  in  the  representation  of  Uq.  For  m  >  7,  we  intro¬ 
duce  further  smootheners  and  proceed  in  the  same  manner.  The  proof  is  completed.  <1 


3  Finite  element  solution 

In  this  section  we  analyze  the  /i-version  of  the  Galerkin  finite  element  method  with 
approximation  order  p.  We  first  identify  the  finite  dimensional  approximation  space 
14  C  prove  two  properties  of  the  best  approximation  in  .ff^seminorm. 

After  an  outline  of  the  finite  element  solution  procedure  (subsection  3.2)  we  turn  to 
the  analysis  of  a  h-p  Galerkin  finite  element  method  for  the  numerical  solution  of  the 
model  problem.  In  subsection  3.3.  we  prove  a  result  on  the  phase  difference  between 
the  exact  and  the  numerical  solutions.  In  subsection  3.4.  we  give  several  stability 
estimates  for  the  numerical  solution  uje-  We  prove  a  discrete  inf-sup-condition  for 
general  p  and  conclude  stability  of  uje  w.r.  to  L^-data,  measured  vn.  L  -  ox  H  -norm. 
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respectively.  We  then  proceed  to  the  proof  of  a  dual  estimate  in  still  lower  norms,  given 
for  specific  ’bubble’  data  as  it  is  encountered  in  the  error  estimation  for  p  >  1.  AH 
stability  constants  are  explicitely  computed  w.r.  to  the  wavenumber  k.  The  analysis 
of  the  finite  element  method  is  concluded  with  error  estimates  (subsection  3.5.). 

3.1  The  approximation  space  Vh 

Assume  that  the  solution  domain  fl  has  been  uniformly  divided  into  n  disjunct  intervals 
A,-  =  (x,_i,x,)  called  finite  elements.  Let  h  =  Xi  —  a;,_i  (stepwidth).  The  set  of  nodal 
points 

=  {0  =  . . .  ,a:„  =  1}  (3-1) 

will  be  called  finite  element  mesh.  We  will  consider  mesh  functions  defined  on  Xh 
and  refer  to  them  by  subscript  ft;  for  the  nodal  value  of  a  mesh  function  Uh  in  a  node 

Xi  G  Xh  we  will  write  shortly  u,-  :=  u/i(i,). 

Let  p  >  0  be  integer  and  A  a  finite  element.  We  denote  by  iS'*’(A)  the  linear  space 
of  all  polynomials  with  domain  of  definition  A  and  degree  <  p.  For  given  mesh,  we 
define  the  space  of  piecewise  polynomial,  continuous  functions 

Vh  :=  5^0, 1)  :=  {s  €  Hlin),  s{x)  €  5p(A.),  i  =  1 . . .  n}  . 

Thus  by  definition  14  C  If  not  stated  otherwise,  we  will  assume  that  14  is 

equipped  with  the  -seminorm. 

Within  the  elements  we  introduce  the  local  coordinate  (  by  linear  mapping  A,-  -> 
/  =  (—1, 1).  The  polynomials  in  S^{I)  are  then  written  as  linear  combinations  of  the 
nodal  shape  functions 

jv.(0  =  .  MO  =  ^  ;  -1  < « <  1 

and  (if  p  >  1)  the  internal  shape  functions 

N,{0  =  4>i-i{0  ;  /  =  3,4,...,p-bl 

where  <f)i  is  written  in  terms  of  the  Legendre  polynomials  Pj 

MO  = 

(see  [SB,  pp.  38/39]).  The  internal  shape  functions  vanish  at  the  element  boudaries, 
forming  the  subspace  S^{I)  =  span  {Ns,  N4, . . . ,  Ap+i}  C  S^{I). 
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Let  Aj  =  be  a  finite  element  of  length  h.  Denote  by  (-jOa  inner 

product  on  A  and  by  ||  •  ||a  the  induced  local  L^.norm.  Similarly  we  use  the  notation 
II  ■  11/  on  /  =  (-1, 1).  On  5P(/),  an  inverse  inequality  (in  p)  is  given  by  the  well  known 
Markov  Theorem^  stating  for  y  €  S^{I)  the  inequality 

lly'll/ <  C'in„(p)p^||p||/  (3.2) 

with  CinM  =  (p+  1)V(P'*V2)  [Be].  Hence  CiUp)  <  4\/2  and  a„v(p)  l/\/2  for 
p  — >  oo. 

AVe  now  prove  that  for  any  u  €  (D)  one  can  find  a  piecewise  polynomial, 

nodally  exact  (on  Xh)  function  s  €  14  s.t. 

•  5  is  an  optimal  (in  h  and  p)  approximation  of  u  in  the  H  -seminorm  and 

•  the  integrals  of  s  are  nodally  exact,  quasioptimal  approximations  of  the  integrals 
of  u. 

Theorem  3.1  (Approximation  in  Y/,):  Let  l,p  be  integers  with  1  <  I  <  p  and  let 
u  e  1).  There  exists  an  s  e  Vh  =  Sl{0, 1)  s.t. 

1.  (nodally  exact  approximation) 

yxi  e  Xh  :  5('")(a:,)  =  u(’"^(x,) ,  m  =  -p  -f- 1, . . . ,  0  (3.3) 

2.  (order  of  approximation) 

(L  \  /-m+1 

— j  |u|/+i,  m  = -p-f  1,...,1  (3.4) 


where  Ca  satisfies: 

1.  (7a(— 1)  =  1  (formal  definition) 

2.  Ca{0)  =  1, 

3.  Ca  decreases  for  0  <  /  <  y/p, 

4.  Ca  increases  for  I  >  y/p  and 


^This  inequality  is  usually  given  in  the  I°°-norm  (cf.,  egs.,  [N,  p.  124]).  We  use  an  L^-variant  [Be]. 
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Ca{p)  =  (0  M 
is  the  maximum  of  Ca{l)  over  I  €  {0, 1, . . .  ,p}. 


Remark  5:  With  respect  to  the  stepwidth  h,  the  estimate  (3.4)  is  the  standard  ap¬ 
proximation  result  of  finite  element  theory  (see,  e.g.,  [Sch,  pp.  46—49]). 

Remark  6:  For  /  =  0, 1  the  statements  are  proven  in  [BKSj.  The  following  argument 
is  a  generalization  of  this  proof. 

Proof  :  We  start  on  the  local  level.  Let  A.-  be  a  finite  element  and  let  I  =  (-1, 1). 
We  write  u'(0  €  H‘{I)  C  L^{I)  as 

OO 

u'(()  = 

where  P,(0  Legendre  polynomials  of  order  i  and  equality  is  understood  in  the 

L^-sense.  Set 

AO 

t=0 

and  define  the  integrals  {i  =  0, 1,2, . . .): 


=  «^“'^(1)  —  J  ‘■'■^^(r)dr 

(3.6) 

=  u(“‘)(l)  -  s("’+^^(T)dr. 

(3.7) 

We  will  now  prove  that  (3.3)  holds.  Let  first  i  —  0,  then  from  eqs  (3.6,  3.7)  we  have 
trivially  u(l)  =  5(1).  Further,  by  definition, 

u(0  =  u(l)  -  u'{T)dT  =  u{l)  -  J2aj  Pj{T)dT  =  u(l)  -  2ao 

=  u(l)  -  r  s'{t)dt  =  u(l)  I  Pj{T)dT  =  s(-l). 

•'-l  j=0 
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Now  we  integrate  using  the  well  known  relation 


to  obtain 


K(f)  =  u(l)  +  a„(/>.(f)  -  P,(f)) 


Integrating  once  more  (we  write  U 

/I  2ai 

^  u{T)dr  =  C/(l)  -  2u(l)  -  2ao  -  — 


Obviously,  the  same  result  is  obtained  from  the  integration  of  the  polynomial  s{()  since 
only  the  coeffient  of  Po  influences  the  result  of  integration  over  the  whole  interval  1. 
By  similar  argument  we  conclude  that  integration  of  the  polynomial  s  on  the  one  and 
the  function  u  on  the  other  hand  leads  to  the  same  result  exactly  p  —  1  times.  Indeed, 
by  replacing  repeatedly  Pi{T)  =  (P/_i('r)  —  1)’  that  with  the  zth 

succesive  integration  of  u(^)  or  5(^)  the  coefficient  ctj  enters  the  set  of  coefficients  mul¬ 
tiplying  Pq*  Since  the  norms  of  and  depend  only  on  the  coefficient  of  Po,  both 
norms  are  equal  until  Po  is  multiplied  by  flp-i,  i.e.  in  general  HO 

and  Thus  nodal  exactness,  eq  (3.3),  is  proved  on  an  arbitrary 

element  and  hence  it  holds  globally. 


Let  us  now  prepare  the  proof  of  estimate  (3.4).  With  above  definitions,  the  error 
of  approximation  is 

e'(0  :=«'«) = 


l=P 

and  from  the  orthogonality  property  of  the  Legendre  polynomials  we  have 

oo  c\ 


(3.8) 


It  can  be  proven  (see  [BKS],  chapter  3)  that  s'  is  the  best  ^^.approximation  to  v!  on 
I  and  the  estimate 

||u'  —  s'll  <  -  |u|/+i  (3-9) 

holds  for  0  <  /  <  p,  where  the  constant  Ca  has  the  properties  2.  -  5.  given  in  the 
theorem. 
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Integrating  the  error  e',  we  get 


j  00-1  °°  a- 

=K)  =  jf  {«'(()-*'«))*  =  go- j(  mdt  = 


After  reordering, 


e(0=  E  +  + 

i=p+l  ^ 


h  — 

'  “  2z  +  3  2i  -  1 


and  the  norm  is 


flp+i  2 


1,2  ^  2  ..  .  K  2  gp+i  2 

^  "  .=^1  2i  +  1  (2p  +  1)2  2p  -  1  ^  (2p  +  3)2  2p  +  1 


We  apply  the  relation  (a  -  6)2  <  2a2  +  2^2  to  obtain  (for  i  >  p  +  1) 


A2  <  2ati  ,  2a?+i 

•  -  {2i  -  1)2  {2i  +  3)2 


and  thus 


■?-  «  2  ^  1  ^  X  J-  V  .P  ^ 

,  S 2i  +  1  -  2,,=  5  2i  +  1  +  ,  •  2i  +  1 


(3.10) 


holds.  Taking  now  into  account  the  second  and  third  member  in  the  r.h.s.  of  eq  (3.10) 
we  get 


1  °°  2 
2  <  lya2_£_ 

-P'^k  ‘2i  +  l 


and  hence 


ell  <  ille'll. 


From  eq  (3.9)  it  then  follows  that 


«ll  ^ 


(3.11) 


(3.12) 


holds  for  1  <  /  <  p. 
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Let  us  conclude  the  local  analysis  showing  an  orthogonality  property  for  e.  Since 
s'  is  the  L^-projection  of  u'  on 

=  0  (3-13) 

holds  for  m  =  0, 1, ...  ,p  -  1.  We  claim  that  e(0  =  -  s'{t))dt  is  orthogonal  to 

5P-2(/). 

Indeed,  for  /  >  0  we  compute 

=  jj(y;(t)-s'(t))f_Cdi)dt 
=  f\u'(t)-s‘(t))(r»  +  \)dt. 

m  +  1  J-\ 

which,  together  with  eq  (3.13),  proves  that  e±S'^“^(/). 


This  completes  the  local  analysis.  By  back-transform  7  A  and  summation  over 
the  elements  we  conclude  eq  (3.4)  for  H^-  and  T/^-norm  ,  i.e.  cases  m  =  1,0  in  eq 
(3.4). 

It  remains  to  prove  eq  (3.4)  for  dual  norms.  We  apply  a  standard  argument  [Sch]. 
By  definition,  for  m  >  1, 

II  II 

e  =  sup  -j-j— . 

|u|m 

Let  P^v  e  5^(0)  be  the  L^.projection  of  u  €  on  5';^(n).  Then  by  orthogonality, 
as  proven  in  step  1, 

(e,v-F’"-h) 

||e||.„=  sup  - rj - 

holds  for  1  <  m  <  p  -  1.  Applying  Schwartz  inequality  and  eq  (3.12)  we  conclude  for 
1  <  /  <  p  and  1  <  m  <  p  —  1,  the  estimate: 


where 


/  1  \  1+1  /  /i  \ 

s||-.  <  C.(()  (-)  I4+.  C.W  (2^)  |^<C(/,m) 

C(l,m)  =  C'a(/)Ca(m)  <  Q)  (’rp)"^/^ 


/+m+l 


This  completes  the  proof  of  Theorem  3.1.  <1 
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3.2  The  finite  element  method 

Let  K  C  H}  (il)  =:  V  be  the  approximation  space  introduced  in  the  previous  subsec¬ 
tion.  As  usual,  function  u  £  Vh  \s  called  the  finite  element  solution  of  the  VP  (1.4) 

if 

Vu  €  V/i  :  B{ufe,v)  =  if,v).  (3.14) 

The  approximation  space  V,  =  5^(0  can  be  written  as  a  direct  sum  of  two  sub¬ 
spaces,  namely, 

where  SI  is  the  space  of  continuous,  piecewise  linear  functions  and 

i=i 

with 

^^(Aj)  =  span{A^3,...,Ar^+i}- 

Here,  are  the  finite  elements,  hence  5p(  Aj)  are  local  spaces  of  "bubble”  polynomials. 

Writing  now  «/e  =  Uh  +  Up,v  =  Vh  +  Vp,  where  Uh,Vh  €  31(0,)  and  Up,Vp  € 
we  have  from  eq  (3.14) 

yvh£Sl{Q):  B{ufe,Vh)  =  {f,vh).  (3-15) 


ctnci  /o  1 

Vup  €  SliCl)  :  5(u/e,  Vp)  =  (/,  Vp).  (3.16) 

From  eq  (3.15)  we  conclude 

Vtij.  e  Slin)  :  B(u/„ui)  =  (/  +  (3.17) 

On  the  other  hand,  eq  (3.16),  inforce  of  S;(tl)  =  fflrS;(A),  decouples  into  n  indepen- 
dent  local  systems 

yweS’^iA):  Ba{up,w)  =  if  +  k\h,w)A,  (3-18) 

where  Ha  is  the  restriction  of  the  form  B  to  the  finite  element  A.  By  formally  solving 
these  equations,  we  express  Up  in  terms  of  /  and  Uh-  The  result  can  be  inserted  into 

(3.17),  giving  rise  to  (3.19) 

Discretising  this  equation  by  standard  approach,  we  obtain  the  linear  system 

[Lfc]{«/i}  =  {r/t}, 


(3.20) 
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where  [Lh]  is  a  (n  x  n)  matrix,  usually  called  the  condensed  stiffness  matrix,  and 
{u/i}  =  Ufe\x^  is  the  vector  of  nodal  values  of  the  finite  element  solution  on  the  mesh 
Xh-  The  piecewise  linear  part  of  u/e  is  determined  by  eq  (3.20),  provided  [Lh]  is 
nonsingular.  The  internal  part  Up  can  then  be  found  locally  by  eqs.  (3.18). 

For  the  sake  of  further  analysis,  we  outline  below  the  details  of  the  solution  proce¬ 
dure; 


Step  1  (Local  approximation  and  static  condensation):  On  any  element  Aj,  the  trial 
function  u  and  the  test  function  v  are  written  as  scalar  products  of  shape  functions 
|A^/,  iVj ,  •  •  •  >  vectors  of  unknown  coefficients  {a^}  =  ^a{,a{, . . . ,  flp+i} 

and  {IP}  =  {b{,bi,...,b’p+iY ,  respectively.  Identifying  a{  =  u{xj-i),a^2  =  “(^i) 
b{  =  vj.rM  =  Vj,  we  have  Uh\^j  =  o,{nI  -f  a{Nl.  The  condition  that  uje  be  the 
solution  of  the  VP  (1.4)  for  all  u  G  14  leads  locally  (i.e.  on  Ay)  to 

(3.21) 

with  the  (p  +  1)  X  (p  +  1)  square  matrix  (/,m  =  1, . .  .,p+  1) 


[Bi]  =\{f  Nj(x)'  Ni(x)'dz  -e  f  Ni(x)  Ni,{x)dx  -  ikN;(\)Ni(l)]  (3.22) 

L  ''^3  ''^3 


and  the  right  hand  side 


}  =  {(/{*),  Nl{x)U, ,  1  =  1, . . . ,  p  + 1}’’ . 


Now,  decomposing 


(B^]  = 


[b;.]  [Bfr] 

K]  [b4] 


(3.23) 


where  is  the  left  upper  2  x  2-submatrix  of  [B^],  and  assuming  for  the  moment 

that  \bL\  is  nonsingular,  we  define 

\CB’]  =  [bJi]  -  [BJs]  [b4]'’  [B5.]  .  (3.24) 

Then,  by  local  variation  of  {^4, . . .  we  find  -  cf.  eq  (3.18)  - 

{Ui-,  u,}  {CB‘\  {x}  = 
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where 


V } = { 1 ; 


(3.26) 


On  uniform  mesh,  the  local  matrices  [C’5'^]  are  identical  on  all  elements  and  can  be 
written  in  the  form 

,CB]  -  [  1  (3.27) 

where  Sp{K)  and  Tp{K)  are  rational  polynomial  functions  of  iif  =  khf2. 

Remark  7:  The  analogy  to  the  /i-version  with  piecewise  linear  approximation  is  given 
by  the  following  consideration.  On  Aj,  the  homogeneous  /i-p- finite  element  solution  is 
written  as 

uye(x)  =  U/e(x,-l)iVf(x)  +  Ufe{xj)NUx)  (3.28) 

where  the  condensed  shape  functions  NJ  €  S^iAj)  are  local  variational  solutions  to 
the  homogeneous  Helmholtz  equation  with  appropriate  boundary  conditions. 

The  local  stiffness  matrix  is  written  as^ 


[CB^] 


B{NIN!)  BiNlN^,) 
BiN^Nl^)  B{N^,N!) 


which  is  again  the  formal  analogon  to  the  /i- version  (cf.  part  I). 


Step  2  (Global  assembling  and  solution  for  Uh):  Enforcing  continuity  of  the  test  func¬ 
tions  in  the  nodal  points  of  Xh  we  obtain  the  set  of  linear  equations  (3.20).  The  discrete 
operator  T/,  is  an  n  x  n  tridiagonal  matrix: 

r  2Sp{K)  Tp{K) 

Tp{K)  2Sp{K)  Tp{K) 

[L.]=  ••.  ,  (3.29) 

Tp{K)  2Sp{K)  Tp{K) 

Tp{K)  Sp{K)-iK 


^see  next  subsection  for  the  proof 
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the  global  right  hand  side  vector  is 


(3.30) 


and  K  =  it/i/2  is  a  measure  for  the  number  of  elements  per  wavelength  (see  part  I, 
remark  8). 


Step  3  (Local  ’’decondensation”):  The  discrete  equivalent  of  eqs  (3.18) 


(3.31) 


can  be  inverted  -  provided  B^2  regular  -  to  determine  locally  the  bubble  modes  of 
the  finite  element  solution. 


3.3  Discrete  Green’s  function  and  discrete  wavenumber 

The  nodal  values  of  the  finite  element  solution  on  Xh  are  found  from  the  tridiagonal 
linear  system  (3.20).  On  uniform  mesh,  this  system  consists  of  formally  identical 
difference  stencils 

Tj,{K)uh{xj.^)  +  2Sp{K)uh{xj)  +  Tj,{K)uhixj+^). 

Equating  these  stencils  to  zero  we  find  the  homogeneous  solutions 

yhi  =  exp{ik'xh),  yh2  =  ex^{-ik'xh)  (3.32) 

where  the  parameter  k'  is  determined  as  a  function  of  fc,  h  and  p  by 

cos{k'h)  =  (3.33) 

The  discrete  wavenumber  k'  is,  in  general,  different  from  k,  causing  a  phase  difference 
between  the  finite  element  and  exact  solutions. 
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The  discrete  finite  element  solution  {u/j}  can  be  written  as 

^  Gh{xi,  Sj)Rhisj)  (3.34) 

i=i 

with  the  discrete  Green’s  function  (cf.  part  I) 


Gfilk  ,  Xj ,  5j) 


1 

sin  k'h 


sin  k'xi  {A  sin  k'sj  +  cos  k'sj) 
sin  k'sj{A  sin  k'xi  +  cos  k'xi) 


X  <  s 
5  <  X  <  1 


(3.35) 


Note  that  the  Green’s  function  does  not  depend  directly  on  p.  The  dependence 
occurs  only  implicitely  through  the  parameter  k'.  This  means  that  the  estimates  of 
G  that  are  given  in  part  I  for  piecewise  linear  approximation  carry  over  to  higher  p 
without  modification  -  except  for  the  estimation  of  k'  in  terms  of  k,  h  and  p. 


Theorem  3.2  (Phase  difference):  Let  p>l  and  k'  be  the  parameter  in  the  funda¬ 
mental  system  (3.32)  of  the  set  of  linear  equations  (3.20). 


Then,  if  hk  <  1, 


(3.36) 


where  k  is  the  exact  wavenumber,  Ca  is  the  approximation  constant  from  Theorem  3.1. 
and  C  does  not  depend  on  k,  h  and  p. 


Remark  8:  The  phase  difference  between  the  exact  and  the  finite  element  solution 
has  been  extensively  studied  in  [HH]  (for  p  =  1)  and  [TPl]  (for  p  =  1,2,3).  As  a 
conlusion  from  numerical  experiments,  the  statement  of  the  theorem  has  been  induced 
in  w.r.  to  k  and  h  in  [TPl].  We  now  prove  that  with  increasing  p  the  phase  difference  is 

also  going  down  with  a  factor  (f (7rp)-^/2(2p)-2p,  i.e.  the  improvement  with  higher 
approximation  is  still  more  significant  than  it  was  assumed  in  above  named  references. 


We  will  give  the  proof  of  this  theorem  after  the  following  preliminary  discussion: 
First  we  observe  that  any  homogeneous  solution  to  the  Helmholtz  equation  can  be 
written  on  each  inner  element  Aj  C  H  as 

u{x)  =  u^ti{x)  +  u^t2{x)  (3.37) 

where  u{xj—\),u^  :=  u{xf)  are  the  nodal  values  of  ti  on  Xh,  whereas  t\,t2  satisfy 

t"  +  =  0  on  Aj  (3.38) 
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with  inhomogeneous  local  Dirichlet  data 

ti{xj-i)  =  l,  ti{xj)  =  0 

(3.39) 

or 

(3.40) 

t2{xj-i)  =  0,  t2{Xj)  =  1, 

resp. 

By  discrete  evaluation  of  the  VP  (1.4)  in  the  nodal  points  of  Xh  we  find  that,  for 

;■  =  -  1, 

To{K)u{xj+x)  +  2So{K)u{xi)  +  To{K)u{xj.x)  =  0 

holds  with 

To{K)  =  =  B{h,tx) 

2So{K)  =  +  B{t2,t2)- 

The  fundamental  solutions  are 

Zhi  =  exp{ikxh),  Zh2  =  exp{-ikxh) 


and 


cos{kh)  = 


So{K) 

UK) 


holds. 

Second,  writing  the  finite  element  solution  on  A_j  as 

Ufe{x)  =  u],N^{x)  +  u%N^{x) 

where  €  S^{Aj)  are  approximate  solutions  to  eq  (3.38)  with  boundary  condi¬ 

tions  ((3.39)  or  (3.40),  resp.,  we  have 

T,  =  BiNlN^)  =  BiN^,N!)  (3.41) 

2S,  =  B{N!,N!)  +  B{NIN^).  (3.42) 


To  see  this,  let  us  analyze  the  BVP’s  (3.38,  3.39)  and  (3.38,  3.40)  on  /  —  (  1,1). 
After  linear  transformation  Aj  I  we  arrive  at 


t"(0  +  kH(()  =  0. 


(3.43) 
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We  assume  that  K  <a<2T.  The  boundary  conditions  are 

t{-l)  =  1,  t{l)  =  0  (3.44) 

or 

t(-l)  =  0,  <(1)  =  1  (3.45) 

resp. 

To  formulate  an  equivalent  variational  problem,  we  write  the  admissible  functions 
as 

t  =  Ni  +  (j)i  or  t  =  N2  +  (f>2 

resp.,  where  (pi  G  H^(^I^  and  Ni^N2  are  the  linear  shape  functions.  The  objective  is 
then  to  find  G  Hl{I)  such  that 

Vt  G  Hl{I)  :  r)  =  -BK{Ni, r)  (3.46) 

holds  for  i  =  1,2,  resp.  (cf.  [SB,  pp.  16/17]). 

The  real  bilinear  form 

Bk{u,v)  =  {u',v')i  -  K^{u,v)i 

is  symmetric  and  coercive  (cf.  Lemma  2.2).  We  can  find  uniquely  defined  ’one  element 
solutions’  by  solving:  Find  <f>^i  eVh  =  S^il)  s.t. 

VctGH:  BK{<f>l<T)  =  -BKiNi,a)  (3.47) 

holds. 

We  now  show  eqs  (3.41,  3.42).  Writing  (with  the  usual  summation  convention  for 
j  =  3, . . .  ,p  +  1) 

N?iO  =  Nr{0  +  a{Nj{0 

NUO  =  N2{0  +  4Nj{0 

we  find  from  eq  (3.31)  the  vectors  {ui},  {02}  to  be 

{ai}  =  -  [522]''  [B21]  I  J  I ;  {02}  =  -  [522]''  [B21]  I  J  I .  (3.48) 

It  is  now  easy  to  see  that 

IBk(K,n;)  Bk(nini) 

[CB]  = 

[Bk(NIN^)  BK(NlNi) 


(3.49) 
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Indeed, 

Bk{NI,  Nl)  =  Bk{NuN^)  +  N^)  +  BK{a\Nu  a'^Nm) 

=  B;,(7Vi,iVi)-2(|  J  I  [5i2])([522rM52i]|  J  |) 

+  (|  J  [B,,]  ([B32I-' IB21]  I  J  }) 

=  B„[l,  >1  -  I  J  }  l^'^l  [522)’'  [B21I I  J  } 

=  CB[1,1] 

and  so  forth  for  C'2B[1,2],C'5[2,1]  and  C'B[2,2].  In  the  equation  chain  above,  we  have 
repeatedly  used  both  the  symmetry  of  the  form  Bk  and  the  local  stiffness  matrix. 

This  validates  eq  (3.49).  Transforming  back  to  global  coordinates  and  assembling, 
we  obtain  eqs  (3.41,  3.42).  <1 

In  the  proof  of  Theorem  3.2.  we  will  also  use  the  following  result. 

Lemma  3.1  Let  ti,<2  and  be  exact  and  finite  element  solutions  to  eq  (3.43) 

with  Dirichlet  data  (3.44)  or  (3.45),  resp. 

Then  for  f  =  1, 2;  j  =  1, 2 

BK{t<  -  K,ti  -  NJ)  =  BK(Nf,Nf)  -  (3-50) 

holds. 

Proof  : 

Let  us  specify  in  eq  (3.46) 

f  =  l,r  =  :  BK{Ni,<t>t)  =  -BK{<f>\.<t>\) 

i  =  2,T  =  rx-  bk{N2,  ri)  =  -BM,  ri) 

i  =  1,T  =  ^2  ’  BK{N\,<f>°2)  =  —BKi<l>i,<l>2) 


(3.51) 

(3.52) 

(3.53) 
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and,  similarly,  in  eq  (3.47) 


11 

b 

II 

Bk{^iA\)  =  -BK{<i>\A\) 

(3.54) 

. 

II 

II 

Bk{^2i^\)  —  ~Bk{4>^i^\) 

(3.55) 

Bk{N\i^2)  —  ~Bk{4'1i4^2)- 

(3.56) 

Furthermore,  since  C  H], 

Bk{^1i4'i)  —  ~BK{<l>ii<fl) 

(3.57) 

i  =  2,  T  =  :  Bk{N2,  <f>i)  —  —Bk{4>2i  *^i) 

which,  together  with  eqs  (3.54)  and  (3.55)  shows 

(3.58) 

BkW 

A\)  =  BMA) 

(3.59) 

Bk{4'2^ 

K)  =  BkW.'K)- 

(3.60) 

We  now  can  check  the  statement  of  the  lemma  by  direct  computation.  For  i  =  j  -  I 
the  l.h.s.  is 

=  Bk{4'\i  4>\)  -  ^BK{<f>i,  4>\)  +  Bk{<I>\,  <I>i) 

=  BK{<i>i,ri)-BK{<f>ifi) 

by  eq  (3.59).  The  r.h.s.  is 

Bk{N!,  Nf)  -  BKiUJi)  =  Bk{Ni  +  (t>\,  Ni  +  (^?)  -  Bk{Ni  +  (j>l,Ni  +  4>\) 

=  Bk{Ni,  Ni)  +  2Bk{Ni,  (j)i)  +  BK{<f>i,  <f>i) 
-Bk{N^,Ni)  -  2BK{N^,c|>l)  -  BKifi,  (!>{) 

=  BK{<f>{,ri)-BK{(l>l,<f>l) 

where  we  have  used  eqs  (3.54)  and  (3.51). 

This  validates  the  statement  for  i  =  ;  =  1.  The  computation  of  the  remaining  cases 
i  =  1,  j  =  2  and  i  =  2,  j  =  1, 2  is  entirely  analogous.  <1 

We  are  ready  to  give  the  proof  of  Theorem  3.2. 

Proof  :  (Theorem  3.2.)  We  will  show  that,  neglecting  higher  order  terms  of 
(kh\ 

I  cos(*:  h)  -  cos{kh)\  <Ck  h  — j  I  ^  I 


(3.61) 
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holds  with  an  independent  constant  C.  Assuming  this  relation  for  the  moment  we 
continue  with  (Ci(p)  :=  CCl{p)f4) 

\  2p 

Let  k'h  -  kh  =  e,  then 

sm(H  +  f)  -  kh 

Since  by  assumption  (^hk ^  1  we  see  that  6  is  small  and  we  may  neglect  higher 
order  terms  in  the  Taylor  expansion  of  2  sin  Thus  we  obtain 


e  <  hkCi{p) 


and  replacing  s  =  k^h  —  kh  we  conclude  the  statement  of  the  theorem. 
We  now  have  to  prove  eq  (3.61)  which  we  write  in  the  form 

Tj,  To  - 


(3.62) 


Consider  any  internal  element  A  being  mapped  on  the  master  element  I  —  (-1, 1).  On 
I,  the  exact  homogeneous  solution  of  the  VP  (1.4)  is  represented  by 


UexiO  = 


where 


sinA'^  cosK^ 
hiO  —  2  sin  A'  ^  2  cos  K 

sinK^  cosK( 

2 cos  a:' 

are  solutions  of  the  BVP’s  (3.43,  3.44)  or  (3.43,  3.45),  resp.,  with  K  =  kh/2. 
The  finite  element  solution  is  written  on  I  as 


-/.(o  =  n),Ni(t) + u),m) 


Finite  Element  Solution  to  the  Helmholtz  Equation  . 


33 


where  are  approximate  solutions  to  the  BVP’s  (3.43,  3.44)  or  (3.43,  3.45),  as 

discussed  in  the  preliminaries. 

By  continuity  of  Bk, 

\BK{tl  -  Nl.trr,  -  7V,„)|  <  (1+  K‘^)\tl  - 
and  applying  now  eq  (3.9)  we  have 

\BK{tl  -  N,,tm  -  iV,„)|  <  (1  +  \tl\p+l\tm\p+l  (3.63) 

P^p 

for  l,m  =  1,2. 

By  direct  computation, 

(  if  p  is  odd 

=  (3.64) 

[  K'^^\ti\l  if  p  is  even 

with 

ll'if  =  I  +  0(K^y.  |(,|f  =  i  +  O(A-). 

Here  and  in  the  following,  0{K^)  means  an  expression  of  the  form  CiK'^  +  C2K'^  +  . . . 
with  constants  Ci  not  depending  on  h,  k  and  p. 

Also  by  direct  computation. 


BK{tutl)  =  BK{t2,t2)  =  \  +  OiK^) 

(3.65) 

and 

BK{tut2)  =  BK{t2,h)  =  -\  +  0{K% 

(3.66) 

Finally, 

we  recall  from  Lemma  3.1  that 

Bxitl  —  Ni,tm  -  Nm)  =  BKiNijNm)  —  BK{ti,tm)- 

(3.67) 

holds  for  /,  m  =  1,2. 

Returning  to  the  proof  of  eq  (3.62),  we  have 

\S,To-SoT,\  =  \BKiti,tr)BK{Nr,N2)-BK{tut2)BKiNi,Ni)\ 


=  \BKih,tl)  {BK{tl  —  Ni,t2  —  N2)  —  BK{ti,t2))  — 

BKiU^h)  {Bk{1\  —  Ni,ti  —  Ni)  —  BK{t\it\))  I 

=  \BK{U-,ti)BK{t\— Ni,t2  — N2) 

—BK{ti,t2)BK{ti  —  Ni^ti  —  A^i)|  (3.68) 


34 


F.  Ihlenburg  and  /.  Babuska 


Here,  we  applied  eq  (3.67)  to  expand  the  expression  on  the  r.h.s.  Thus 

\SpTo-SoTj,\  <  \BK{tuh)\\BK{ti-Ni,t2-N2)\  + 

and  applying  eqs  (3.63,  3.65,  3.66)  leads  to 

\S,To  -  SoTp\  <  (^  +  0{K^)^  ((1  +  {\ti\p+i\t2\p+i  +  l^ilp+i))  • 

For  odd  p  we  then  have  directly  by  (3.64) 

|5,T.  -  5„T,|  <  (i  +  0(K^))  (d  +  (5  +  0(/<’))) 


Caipf 


,J^2p+2 


(3.69) 


where  we  neglected  terms  of  order  0{K^).  Then  also  (3.62),  and  hence  the  statement, 

^  ^  -  \SpTo-SoTA 

Tp  To  \ToTp\ 

and  it  can  easily  be  seen  that 

ir^Tpl  =  \BK{tut2)\\BK{t^  -  -  N2)  +  BK{tl,t2)\ 

is  bounded  from  below  by  a  constant  independently  from  K,p. 

If  p  is  even  then  we  insert  eqs  (3.65,  3.66)  into  eq  (3.68)  to  ge 

l^pTo-^-oTpl  <  ^\BK{h-Ni,ti  +  h-{Ni+N2))\  + 

0{K^)  {\BK{ti  -  Nuh  -  ^2)1  +  \BK{ti  -  NuU  -  Ni)\) 

For  the  first  term  in  this  equation  we  have 

\BK(U-Ni,ti+h-(Ni  +  N2))\  <  l(*i+*2)lrt-i 

—  p2p  COS  K  J 


<  CCaip) 
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with  C  not  depending  on  K,p.  Again,  terms  of  order  0{K^)  have  been  neglected. 
For  the  second  term,  a  similar  estimate  follows  directly  from  eqs  (3.63,3.64).  Thus  the 
estimate  (3.69)  holds  also  for  even  p,  and  the  statement  follows  by  similar  argument.  < 


3.4  Discrete  stability 

In  Part  I,  when  investigating  the  /i- version  with  p  =  1,  we  proved  the  stability  estimate 

llt^'/ell  <  C\\f\\ 

and  showed  that  the  inf-sup-constant  on  the  discrete  subspace  is  ^1^  =  Ck~^.  A  stan¬ 
dard  corollary  then  yields  the  stability  estimate 

\\u'j^\<ck\\r% 


We  will  now  show  that  both  results  carry  over  to  higher  p. 

Further  we  will  define  in  this  subsection  a  specific  data  subspace  that  we  will  en¬ 
counter  in  the  error  analysis.  In  this  subspace  we  will  prove  stability  w.r.  to  higher 
integrals  of  the  data  -  the  discrete  analogon  to  Theorem  2.2. 

Let  hence  u/e  €  be  the  finite  element  solution  to  the  VP  (1.4)  for  data  /  € 

We  write 

Ufe  =  U/i  +  Up 

where  Uh  is  based  on  the  nodal  shape  functions  and  Up  on  the  internal  ones.  From 
the  definition  of  the  shape  functions  Nj  and  orthogonality  property  of  the  Legendre 
polynomials  it  follows  that 

\W,£  =  ll-'if  +  Hf-  (3.70) 

Furthermore,  Up  satisfies  on  each  element  A  eq  (3.18).  Transforming  A  —+/■*■  = 
(0, 1),  we  get 


Viy  €  ^^(A)  :  ^^(up,  w)  =  /i^(/,  w)  ^  K'^{uh,  w). 

We  now  prove  a  first  stability  lemma  on  |lUp||  for  data  /  € 


(3.71) 


Lemma  3.2  Let,  for  uje  be  the  finite  element  solution  to  the  VP  (l.f)  data 
f  G  L'^(Cl).  Assume  that  hk  <  a  <  tt. 

Then 

h'jeW  <  C\\f\\ 

holds  with  a  constant  C  independent  of  h,k  and  p. 


(3.72) 
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Proof  :  Let  «/e  =  Uh  +  Up  as  defined  above.  We  know  (cf.  part  I,  Lemma  3)  by 
straight  estimation  of  the  discrete  Green’s  function  representation  that  ||u/j||  <  GiU/H, 
where  Ci  does  not  depend  on  h,  k  and  p.  Applying  eq  (2.5)  from  Lemma  2.2  and 
Remark  2  to  eq  (3.71),  we  have 

IKII  <  D  (A’li/ii  +  /r^ll^ill) 

where  D  —  7r/(7r^  —  O!^).  Ba-ck-trsinsforin  to  ^  then  yields 

IKIU  <  Dh  (II/IU  +  *:'I|xi||a)  • 

Summing  up  and  applying  Schwarz  inequality,  we  get 

IKII  <  Dh  (ll/ll  +  P||u/.||)  <  \\f\\h  {D  +  C2k) 
where  we  applied  A;||u;i||  <  C'2ll/ll  with  C2  not  depending  on  h^k^p.  Thus 

llu'^J  <  (Cl  ■^Dh  +  Cikh)  ll/ll, 

and  the  statement  is  readily  obtained  by  neglecting  Dh  and  setting  (7  =  Gi  +  G2Q:. 
The  proof  is  completed.  <d 


We  now  prove  the  inf-sup-condition  for  T4  —  ^h(D)' 

Lemma  3.3  (Discrete  inf-sup-condition):  Let  Vh  =  ^^(fi)  and  let  B  iVhxVh 
(E  be  the  sesquilinear  form  defined  by  eq  (1-4)- 
Then,  if  h  is  such  that  hk  <  a  <  ir, 


\B{u,v)\  C 
mf  sup  I  I  I  I  '■  >  -r 


(3.73) 


where  C  does  not  depend  on  h,  k  and  p. 

Proof  :  The  argument  is  similar  to  the  case  p  =  1.  For  arbitrarily  fixed  u  G  Vh  set 
v:=u  +  z  where  z  £  Vh  is  solution  to  the  VP 

VweVh:  B{w,  z)  =  k’^iw,  u). 
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The  ratio  k/k'  is  bounded  for  kh  <  tt  and  p  >  2  by  Theorem  3.2.  Further,  Zp  solves 
locally,  i.e.  for  all  w  €  ^^(A), 

J5A(zp,iy)  = 

Applying  Lemma  2.2  and  Remark  2, 

IKII<CiKM1  +  II^/.II)- 

Again  by  discrete  Green’s  function  representation,  we  can  show  that  ||5r/i||  <  C'2||u  |1, 
with  C2  not  depending  on  h,k  and  p.  Applying  a  Poincare  inequality  for  ||it||,  we 
conclude 

Ill'll  <  CskWu'W 

and  the  statement  is  concluded,  using  the  particular  choice  of  z  (cf.  part  I).  <d 

It  follows  by  standard  theory  that  |uli  <  Ck\\f^~^m. 

We  collect  both  stability  estimates  in  the  following  proposition. 

Theorem  3.3  (Stability  of  the  FE-solution  I):  Let  f  6  and  let  u/e  €  14  = 

be  the  finite  element  solution  to  the  VP  (l-4)‘ 

Then,  if  h  is  such  that  hk  <  a  <ir,  the  stability  estimates 


I«/.I.  <  C, ll/ll 

(3.74) 

(3.75) 

hold  for  Cl,  C2  not  depending  on  h,  k  and  p. 

Next  we  formulate  a  dual  stability  property  that  we  will  use  to  prove  a  preasymp- 
totic  error  estimate  for  the  p-degree  finite  element  solution.  To  this  end,  we  define  a 
specific  data  subspace. 

Definition  3.1  For  integer  l>0  we  define  a  subspace  F^{Q,)  C  L^{Q.)  by 

^  1  U/.  =  0  i  =  1 . . . /} 


with  Ff{n)  :=  L\n). 
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Observe  that  for  this  space  we  can  show,  by  adding  up  local  Poincare  inequalities, 

l|/(-)||</i||/(-’+')||.  (3.76) 

In  particular,  for  /  €  Fj  and  hk  <  a,  inequality  (3.74)  directly  follows  from  (3.75). 

For  the  discrete  data  obtained  from  f  €  Fj(fi)  by  the  finite  element  procedure,  the 
following  proposition  is  true. 

Lemma  3.4  Consider  the  VP  (1.4)  on  14  =  SK^)  with  data  f  €  Let  Aj  be 

an  arbitrary  finite  element  and  let  be  the  condensed  right  hand  side  vector 

given  by  eq  (3.26).  Assume  further  that  the  stepwidth  h  is  sufficiently  small  so  that 
hk  <  a  <  TT. 


Then 

|f,i  <  (3.77) 

holds  for  even  m  =  0, 2, . . .  <  p  —  1  with 

Cd{p,0)  =  l 

C,(p,  m)  =  C.  +  m  >  2 

where  Ci,C2  do  not  depend  on  h,k  and  p. 

Proof  :  Let  first  m  =  0.  We  omit  in  the  notation  the  element  number  j  and  renumber 
formally  {rj-i,fj}  Eq  (3.26)  then  reads 


[B12]  [B22] 


It  is  straightforward  to  show  that  rj  =  (/, A^i)a  ^2  —  {fiLl2)A  satisfy  |rj|  < 
Define  {y}  :=  [522]“^  {^*3,  •  •  •  This  vector  is  the  discrete  solution  to  the  VP 

VuieF^A):  B^{y,w)  =  {f,w)A 

and  by  Lemma  2.2  and  Remark  2, 
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By  simple  computation,  Hj/'IIa  >  hence 

|{y)l  <  C.A=/1/||a, 

where  C\  =  7r/(7r^  —  a^)  does  not  depend  on  h,  k  and  p. 

Consider  now  the  term  {z}  :=  [5i2][-B22]~^{^3> •  •  •  —  [■®i2]{y}-  The  coef¬ 
ficients  of  the  matrix  [5i2]  are  bij  =  Nj)  with  z  =  1,2  and  j  —  3...p-f  1. 

For  these  i,j  we  can,  using  integration  by  parts,  that  =  0  and  hence 

I-.  =  -k^{Ni,Nj)i,.  Obviously,  the  euclidian  norm  of  the  rows  {bi},i  =  1,2  can 
be  bounded  as 

\{bi}\<C2k^h. 

The  constant  C2  does  not  depend  neither  on  h,  k  nor  on  p  since  the  bandwidth  of  the 
local  mass  matrix  does  not  increase  with  p  for  p  >  3  (cf.  [SB,  p.46]). 

Thus,  with  the  previous  estimate  of  |l{j/}|l,  for  z  =  1,2, 

l^.i  =  IhVjl  ^  |{^«}lAl{y}|A 

<  CiC2k^h^^^fU 

<  Csh^^^fWA  (3-78) 

with  C3  =  CiC2a^-  Together  with  the  observation  for  ri,r2,  the  last  estimate  proves 
the  statement  for  m  =  0. 

Let  next  m  >  2  and  assume  for  convenience  that  A  has  been  mapped  to  /'*'  =  (0, 1). 
Then,  for  j  =  1,2, 

fj  =  h  [' f{e)ipj{9)dd  (3.79) 

Jo 

where  <pi,V=’2  €  S^^{I'^)  or  resp.,  are  the  ’one-element’  solutions  to  the  homo¬ 

geneous  Helmholtz  equation 

u"  +  k'^h\  =  0  (3.80) 

with  the  boundary  conditions 

u(0)  =  l,  u(l)  =  0  (3.81) 

or 

u(0)  =  0,  zz(l)  =  l,  (3.82) 

resp.  The  exact  solutions  and  ^2  of  these  BVP  are 

ti{9)  =  cos  kh9  —  cot  kh  sin  kh9 
t2{9)  =  sin  kh9 1  sin  kh. 


(3.83) 

(3.84) 
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Integration  by  parts  in  eq  (3.79)  leads  to 

I?, I  =  I  /' 

JO 

no  boundary  terms  occur  for  m  <  p  —  1  since  /  €  Consequently  we  have  for 

For  the  estimation  of  |l<^j”*^||/+>  we  define  Xj  €  or  resp.,  by 

Xi(0)  :=  ri(0)  -0ri(l) 


X2(^)  :=  '^2(^)  +  ^(1  “  '’^2(1)) 

where  ti,T2  €  are  the  Taylor  polynomials  of  order  p  in  =  0  for  ti{6)  and 

t2{0),  resp.  Now  trivially 

As  to  the  second  member  on  the  r.h.s.,  it  can  be  shown  by  direct  computation  that  for 
even  m  >  2 

llxf^ll/^  <  C^ihkr.  (3.85) 

holds  with  a  constant  Ci  not  depending  on  h,k  and  p.  Turn  to  the  estimation  of 
-  Xj"'^l|/+-  For  m  =  1,  we  have 

|(Pi  -  Xili  ^  l^j  ~  V’ili  +  l^i  “  Xili  ^  C'2|tj  -  Xj|i  (3.86) 

by  Cea’s  lemma.  ^  By  construction  of  Xj, 

K,  -  »|.  =  K,  -  r,l  +  O  (<^)  <  Ca  (^^)  ,  (3.87) 

where  C3  does  not  depend  on  h,k  and  p. 

For  m  >  2,  we  use  repeatedly  the  inverse  inequality  (3.2)  to  relate 

ii^r*  -  xr’ii/*  s 

■^To  be  precise,  we  write  as  sums  of  linear  and  bubble  functions  and  apply  the  statement  of 
Cea’s  lemma  on  an  appropriate  subspace  Vh  C  Ho{I^)- 
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Inserting  this  estimate  into  (3.86)  vs.  (3.87),  we  conclude  that 

llv’S  -  xj  'II/*  <  C,  (“) 

where  Ca  =  does  not  depend  on  h,k  and  p. 

liar'll, + <  c,(hkr + 


and,  finally 


<  h{hhr  (c. + 


The  statement  now  follows  by  back-transform  /+  -»  A.  The  proof  is  completed.  <1 

Remark  9:  For  odd  m  this  statement  holds  true  with  the  additional  assumption  that 
hk  is  bounded  from  below,  i.e.  0  <  ^  <  hk. 

Let  us  show  this  for  m  =  1  and  j  =  2.  We  have 

I|X'21I<IKII  +  111-^2(1)11 

By  construction,  the  second  member  is  O  Also  r^iO)  =  t2{0)+O  )•  Hence, 

neglecting  terms  of  higher  order  (note  that  by  assumption  p  >  m), 

kh 

llx'J  <  ll«'/ll  =  ^ll“=“«ll'* 


Furthermore,  \ip2  —  X2I1  =  ^  Hence  ||y>2ll/+  —  ^2  and 

|r2|  <  C2h\\f~'^^\i+ 

after  back-transform  7+  ^  A.  Multiplying  and  dividing  now  by  we  get 

IF2I  <  ^2^^||/”^^ll/+ 

<  Ch^^^k\\f(-'^% 

with  C  =  C2y9-L  This  shows  our  case;  the  argument  for  j  =  1  or  higher  m  is  similar. 
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Corollary  3.1  For  the  norm  of  the  discrete  right  hand  side  Rh,  the  estimate 

(3-88) 

holds  for  m  =  0 . . .  <  p  —  1. 

Proof  :  By  definition  (cf.  part  I), 

i=l 

where  the  coefficients  Ri  are  given  by  eq  (3.30).  By  Lemma  3.4, 

/  n  \i/2 

and  the  statement  readily  follows.  <1 

We  now  prove  a  proposition  on  dual  stability  if  the  data  f  G  F^  (fl). 

Theorem  3.4  (Stability  of  the  FE-solution  II):  Let  uje  €  14  =  ‘S'^(^) 
solution  to  the  VP  (l.f)  with  data  f  G  F^{n),  where  m  is  even,  m<p-l.  Assume 
that  0  <  kh  <  a  <  TT . 

Then  ,  ^ 

(3.89) 

where  Cd  is  the  discrete  dual  stability  constant  (cl.  Lemma  S-i)  and  Ci  does  not  depend 
on  h,k  and  p. 

Proof  :  In  light  of  Theorem  3.3.  we  only  need  to  prove  the  statement  for  m  ^  2. 
As  before,  we  write  Ufe  =  Uh  +  Up.  Then 

ll^/ell  ^  +  C'i(||/^“^^ll  +  k^\\u[  ^^11) 

by  Lemma  3.2.  It  can  be  shown,  using  the  Green’s  function  representation  of  Uh  that 

ll^lll  ^11  — 

where  C2  does  not  depend  on  k  and  p»  On  the  other  hand, 

||H||  <  G(p,m)AA:”i|/(-”>||, 

from  Lemma  3.4,  and  the  statement  follows.  The  proof  is  completed  <1 

Remark  10:  The  stability  theorem  holds  for  odd  m  with  the  additional  assumption 
that  0  <  ^  <  hk  -  d.  the  previous  remark. 
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3.5  Error  estimates  for  the  finite  element  solution 

Throughout  this  subsection,  we  denote  by  I  the  regularity  of  the  exact  solution  w,  i.e. 

we  assume  that  u  E  :=  fl  V  with  V  = 

In  part  I  we  proved  that  the  finite  element  solution  is  asymptotically  quasioptimal. 
The  same  result  holds  for  higher  approximation.  However,  the  range  of  asymptotic 
behaviour,  as  taken  w.r.  to  the  meshsize  grows  with  p. 

Theorem  3.5  (Asymptotic  estimate):  Let  I  >  l,p  >  1.  Let  further  u  e  and 

Ufe  ^Vh  =  Sl{0.)  be  the  exact  and  finite  element  solutions  to  the  VP  (l.f),  resp. 
Then,  if  k'^hfp  is  sufficiently  small,  the  quasioptimal  estimate 

|u  -  Ufe\i  <  C  inf  |u  —  uli  (3.90) 

holds  with 


provided  the  denominator  of  C  is  positive. 

Proof  :  Denote  e:=u-Ufe  and  let  ^  be  the  solution  to  the  VP 

\/weV:  B{w,z)  =  {w,e).  (3.92) 

This  problem  has  a  unique  solution  2  €  H^{n).  We  can  show  (part  I,  Theorem  3)  that 

\\ef  <  2  (I2  -  u;li|e|i  +  k^z  -  i«I|||e|l) 

holds  for  all  w  e  V.  In  particular,  for  w  =  s,  where  s  is  the  approximation  of  2  as 
constructed  in  Theorem  3.1.,  we  have 

||e|P<2(A|.b|e|.  +  i'(Aj'‘|,l3||el|). 

We  now  apply  the  stability  estimates  (part  I,  Lemma  1,  and  eq  (2.24),  Theorem  2.1) 
to  get 

l^ls  <  (1 +4fc)||e|li  <  ^{1 +4fc)lc|i 


-  ifY  (1 + (“)') 


(3.91) 
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by  Poincare  inequality  since  e  € 
Hence 


iiciicc  2  I — \ 

WP<2(^(l  +  <^)  +  *"(l+4*^)  (^)  ||ell|e|,; 

dividing  now  by  ||e||  and  neglecting  terms  where  h  is  of  order  higher  than  k,  we  arrive 
at  the  intermediary  result 


3  f  kh 


(3.93) 


In  the  next  step  we  use  H-orthogonality  of  e  to  show  (cf.  again  part  I,  Theorem  3) 
that 

l\e\l<6k^e\\^  +  A\u-v\l  +  k^u-vf 

Zi 

holds  for  all  v  €  Vh-  We  choose  v  =  s  from  Theorem  3.2,  then 

h-4  < 

and  applying  now  eq  (3.93)  we  obtain 

i|e|;  -  6P  (y )'  (1  +  4 


hence 


where 


|eli  <  C\u  -  s\i 

-  6^2  (f )"  (1  +  vf  (")") 


This  completes  the  proof.  <1 


Remark  11:  For  piecewise  linear  approximation  we  proved  (part  I,  Corollary  2) 


\u-ufe\  <  Cl  inf  |u-u|i 
veVh 


Ci  = 


2(l  +  (ir)' 

(i-6CfW(l  +  fc)2)' 
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where  2 

Obviously  C2  >  2,  hence  k^h^  <  ill  is  necessary  for  well- definiteness  of  Ci.  A  similar 
computation  yields  k^h^  <  as  a  necessary  condition  for  well-definiteness  of  C  in  the 
theorem  above. 

Remark  12:  In  the  form  given  in  Theorem  3.5,  the  quasioptimal  estimate  holds 
independently  on  the  regularity  of  the  solution  u.  The  order  of  convergence,  in  terms 
of  is  obtained  by  introcucing  the  approximation  property  of  the  subspace  Vh  from 
Theorem  3.1.  For  given  p,  the  maximal  order  of  convergence 

\u  -  u,.|.  <  C  (eftf  (Ej"  |«|,+. 

is  achieved  when  the  regularity  of  u  is  /  >  p. 

We  now  proceed  to  error  estimates  in  the  preasymptotic  range  (i.e.  without  restric¬ 
tions  on  k‘^h).  Let  us  first  relate  the  finite  element  solution  to  best  approximations  in 
Vh  =  Sh(fl),  as  constructed  in  Theorem  3.1. 

Lemma  3.5  Let,  for  p  >l,u  and  uje  be  the  exact  and  finite  element  solutions  to  the 
VP  (1.4),  resp.,  and  let  s  €  Vh  be  a  nodally  exact  quasioptimal  approximation  to  u  in 
the  sense  of  Theorem  3.1. 

Then  z  :=  uje  —  s  is  the  finite  element  solution  to  the  VP  (1.4)  k^{u  —  s). 

Proof  :  Trivially  z  =  Ufe  -  u  +  u  -  s,  and  hy  R-orthogonality  of  u  -  uje  to  14, 
B{z,v)  =  B{u  -  s,v)  holds  for  all  v  €  14-  The  boundary  term  in  B{u  -  s,u)  vanishes 
due  to  =  s\x„.  Using  local  exactness  of  the  integrals  of  s  we  show,  repeatedly 
integrating  by  parts,  that  also  the  term  ((u  —  s)^,  vanishes.  Thus,  for  all  v  €  14> 

B{u  -s,v)  =  -k‘^{u  -  s,  u) 

which  completes  the  proof.  <I 

It  is  now  straightforward  to  show  a  first  error  estimate. 

Theorem  3.6  (Preasymptotic  estimate  I):  Let,  forl<l<p,u£V‘  and  uje  €  14 
be  the  solution  and  the  finite  element  solution  to  the  VP  (I.4),  respectively.  Assume 
that  hk  <  a  <  %. 
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Then  for  e  :=  u  —  u/e 

|e|i  <  Ca{l)  ^1  +  Cik  j  |u|/+i  (3.94) 

holds,  where  C\  does  not  depend  on  h,  k  and  p,  whereas  Ca{l)  is  the  approximation 
constant  (Theorem  3.1),  being  at  most  of  order  . 

Proof  :  Lei  z  =  Ufe-u  +  u-s  BiS  above.  By  Theorem  3.3  and  the  previous  lemma, 
l^li  <  Ck^u  -  s||,  hence 

|e|i  =  |2  +  u  -  s|i  <  Ck‘^\\u  -  s|l  +  |u  -  sji.  (3.95) 

To  complete  the  proof,  we  insert  the  appropriate  results  from  the  approximation  the¬ 
orem.  <1 

Remark  13:  Note  that,  if  k'^h/2p  is  bounded,  the  error  estimate  is  equivalent  to  the 
asymptotic  quasioptimal  estimate  given  in  Theorem  3.5.  This,  again,  is  an  analogy  to 
the  /i-version  with  p  =  1  (part  I,  section  3.6). 

The  estimate  of  the  previous  lemma  can  be  generalized  for  p  >  2,  employing  dual 
stability  estimates  for  the  data  k^{u  —  s)  € 

Theorem  3.7  (Preasymptotic  estimate  II):  Lef  1  <  /  <  p  and  0  <  m  <  p,  m 
even,  with  p  >  2.  Let  u  ^  V‘  be  the  solution  to  the  VP  (1.4)  with  data  f  € 
and  let  uje  €  Vhbe  the  finite  element  solution  to  this  problem.  Assume  further  that  the 
stepwidth  h  is  such  that  hk  <  o:  <  tt. 

Then 

\e\i<Ca{l)  l  +  +kCd{p,m)Ca{m)(^^'^  (^) 

holds  with  Cl  not  depending  on  k,h,p. 

Proof  :  Let  s  e  14  be  a  nodally  exact,  optimal  approximation  of  u  in  the  sense 
of  Theorem  3.1  and  define,  as  before,  z  :=  u/e  —  Uj.  We  know  that  z  solves  B{z,v)  = 
-k'^{u  -  s,v)  for  all  v  G  14.  The  data  of  this  problem  is  in  the  space  Ff-'^,  hence  by 
Theorem  3.4. 

I^li  <  it"  (Cj(p.m)t”||(u  -  s)(-">||  +  Ci||(u  -  s)<-'>||) 
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holds  for  m  <  p  —  1.  Inserting  the  from  Theorem  3.1 

(1  \  /+Tn+1 

we  conclude  the  statement.  <3 

Remark  14:  With  the  additional  assumption  0  <  /3  <  kh,  the  statement  holds  a,lso 
for  odd  m  -  cf.  remarks  9,10.  This  assumption  is  consistent  with  the  error  estimation 
in  the  preasymptotic  range  and  with  computational  application  where  the  magnitude 
of  hk  is,  for  medium  and  high  k,  bounded  from  below  by  practical  considerations. 

Remark  15:  We  obtain  estimate  I  from  estimate  II  by  setting  m  =  0,  hence  II  gener¬ 
alizes  I. 

Let  us  specify  the  estimate  II  for  a  certain  type  of  solutions,  namely,  those  oscillating 
with  frequency  k.  These  solutions  are  of  practical  importance  in  physical  applications 
in  wave  propagation  and  wave  scattering;  they  are,  among  others,  produced  by  Dirac 
data  (point  sources). 

Thus,  having  in  mind  to  specify  solutions  that  essentially  behave  like  exp{ikx),  we 
define: 


Definition  3.2  Let,  for  I  >  1,  u  e  V‘  be  a  solution  to  the  VP  (l.f)-  ca//  u  an 
oscillating  solution  if 

l«|/+i  <  Dk^\u\^.  (3.97) 

holds  with  a  constant  D  not  depending  on  k. 

With  this  definition,  we  directly  have  the  following  corollary. 

Corollary  3.2  (Error  estimate  for  oscillating  solutions):  Let  1  <  I  <  p  and 
p  >  2.  Assume  that  there  is  given  a  data  f  €  H^  such  that  the  solution 

u  €  V'  to  the  VP  (1.4)  is  oscillating.  Assume  further  that  the  stepwidth  h  is  such  that 
hk  <  a  <Tr.  Let  uje  €Vh-  be  the  finite  element  solution  to  the  VP  (1.4)- 


Then  the  relative  error  le|i  :=  [u  —  u/eli/|^|i  is  bounded  by 

'hkV^  /.i\2l  //.LX'+m+l 


e  1  < 


'  ^  (khV 

Cl  +  C2  f  I 

\^pJ 


+  kCs 


:) 


(3.98) 
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where  Ci  =  DCa{l),C2  =  ECaiJ)  ofid  C3  —  DCd{Pti^}Ca{i)Ca{'i^))  with  D,E  not 
depending  on  h,  k  and  p. 

Proof  :  We  introduce  the  definition  of  oscillatory  behaviour  into  the  estimate  II, 
eq  (3.96).  <] 

Let  us  consider  special  cases  of  the  previous  corollary. 

1.  If  /  >  p  we  have  with  0  := 

|e|i  <  e{Ci  +  +  C3k9\ 

This  is  principally  the  estimate  that  was  given  in  the  analysis  of  the  h-version  - 
cf.  Introduction,  eq  (1.5)  -  with  6  —  kh.  Note  also  that,  formally,  the  error  is 
written  as  the  sum  of  best  approximation  error  plus  pollution  term  of  the  order 
or  the  phase  lag.  However,  the  constant  C3  depends  on  p  (see  next  section). 

2.  In  the  case  of  lower  regularity  (/  <  p)  the  pollution  is  for  higher  p  relatively  (i.e. 

compared  to  the  best  approximation  order)  still  smaller  as  in  the  case  of  full 
regularity.  Consider  the  lowest  possible  case  of  Dirac  data.  Then  /  1  and  the 

estimate  is 

|e|i  <  e{Cx  +  020'^)  +  C3k0^^'^ 

In  general,  the  constant  C3  depends  on  m.  Note  that,  for  fixed  approximation 
order  p,  one  is  free  to  choose  m  in  the  range  of  0, ...  ,p  —  1.  This  can  be  used  to 
optimize  the  size  of  the  pollution  term  C3(p,  ?7i)(/iA;/2p)^'*'’”'^^. 

Remark  16:  We  will  show  in  the  numerical  evaluation  that  the  constants  €1,62  are 
sharp.  On  the  other  hand,  the  theoretically  predicted  growth  in  C'3(p)  was  not  ob¬ 
served  in  several  numerical  examples.  It  is  an  open  question  whether  the  dual  estimate 
is  sharp  in  the  pollution  term. 

We  conclude  this  subsection  with  an  error  estimate  in  negative  norms.  We  first 
show  a  mapping  property;  the  proposition  then  readily  follows. 

Lemma  3.6  Let  u  e  V‘  and  u/e  E  Vh  be  the  exact  and  the  finite  element  solution  to 
the  VP  (1.4)-  Then,  fori  <m<p-l, 

||e||-„<DC,(m  +  l)S:^^^  |e|.  (3.99) 

where  Cs  is  the  stability  constant  from  Theorem  2.1  and  D  is  a  constant  not  depending 
on  h  and  k. 
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Proof  :  Since  X  =  H^{n)  is  a  Hilbert  space  there  exists  Vo  E  X  s.t. 

II  II  _  l(e,'>o)l 

"  "  |t>.U  ■ 

Let  ^  U  be  the  solution  of  the  VP  (1.4)  with  data  v^.  Then  by 

Theorem  2.1.  we  have 

\zU2<Cs{m  +  l)n\Vo\\m 

where  Cs  grows  at  most  linearly  with  m.  On  the  full  norm  ||  •  ||m  is  equivalent  to 
the  seminorm  |  •  |m,  hence  there  exists  a  constant  Ci  s.t. 

k|m+2  <  CiC,{m  +  l)k^\vo\m.  (3-100) 

Further,  from  B{z  -  X, «)  =  (^07  e)  for  all  x  €  V/i  we  conclude 

||e||-mbo|m  <  Co(A;)^infJ^-xlile|i 

<  ce  \zU^\e\^ 

where  the  continuity  property  of  the  form  B  and  the  approximation  property  of  14 
have  been  used.  The  statement  now  follows,  inserting  eq.  (3.100).  < 

Theorem  3.8  (Dual  error  estimate):  Let  u  €  and  Ufe  ^  Vh  he  the  exact  and 
the  finite  element  solution  to  the  VP  (1.4).  Assume  1  <  m  <  p  —  1  and  1  ^  I  ^  p. 
Then,  if  hk  <  a  <  t , 

||e||-.  <  C(m,  /)  [C,  +  g j  ||/I!p-i  O-lOl) 

holds,  where  C{m,l)  =  Ca{l)Cs{l)Csim  +  l),  whereas  Ci^C^  do  not  depend  on  h,k  and 
P- 

Proof  :  Combining  Lemma  3.7  with  the  preasymptotic  estimate  I,  as  given  in 
Theorem  3.6,  we  obtain  an  estimate  w.r.  to  luj/+i.  The  statement  is  then  concluded 
from  Theorem  2.1  <1 

Remark  I7j  All  propositions  on  the  finite  element  solution  contain  the  assumption 
hk  <  TT.  Essentially,  we  ensure  herewith  that  the  inversion  of  the  local  stiffness  matrices 
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is  well  defined  {hk  has  to  be  smaller  than  the  minimal  eigenvalue  of  the  condensation). 
However,  the  inversion  is  also  well  defined  if  hk  lies  between  the  first  and  the  second 
eigenvalue  and  so  forth  (to  be  more  precise:  is  bounded  away  both  from  the  first  and 
second  exact  and  the  corresponding  numerical  eigenvalues).  It  is  expected  that,  under 
this  condition,  more  than  one  halfwave  can  be  resolved  by  one  element  of  higher  ap- 
proximation  in  practical  computations. 


4  Numerical  evaluation 

In  this  final  section  we  present  some  results  of  numerical  experiments  that  illustrate 
the  theoretical  results  obtained  in  the  previous  section.  We  solve  the  model  problem 
for  the  constant  data  /  =  1.  The  exact  solution 

u[x)  =  -^  ((1  —  cos  kx  —  sin  fc  sin  kx)  +  i  (1  —  cos  k)  sin  kx)) . 
is  regular  and  oscillating. 

Best  approximation  in  In  Theorem  3.1  we  constructed  a  function  s  G 
that  has  the  best  -approximation  property  to  a  given  function  u  €  /f;+i(fi).  Writing 
s  locally,  after  scaling  to  I  =  (-1,1),  as 

*(«  =  E 
1=1 

we  have,  by  interpolating  property,  =  Uj_i  and  (^2^  =  uj.  Further,  by  orthogonality, 
cP  =  (u'(^'>,iV/).  and  \eba\l  =  \u-  s\l  =  |u|2  -  \s\l  with 


j=l  \  ^  »=3 


The  approximation  property,  as  given  in  Theorem  3.1,  reads  for  m  =  1  and  /  -  p  as 

le-l.  :=  ^  <  C.(p)  (y '  ^  (4.1 


where  Ca  is  the  approximation  constant,  growing  as  (e/2)P/p(l/4). 
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In  Fig.  1  the  error  |e|i  is  plotted  for  k  =  50  and  hk  <  tt.  The  slopes  are  p,  in 
accordance  with  eq  (4.1). 

Consider  next  a  plot  for  the  approximation  error  on  coarse  grid  (/ifc  >  tt  -  Fig.  2). 
We  observe  an  interesting  effect  whenever  the  mesh  locally  hits  the  size  of  a  halfwave 
(two,  three,  . . .)  halfwaves.  To  explain  this  effect,  note  that  for  k  =  mr  (m-integer)  the 
solution  u  reduces  to 


Hence  in  this  case  we  locally  solve  the  problem  of  best  if  ^-approximation  of  uioc  =  cos  x 
on  one,  two,  three,  . . .  halfwaves.  For  even  n,  this  approximation  problem  is  symmet¬ 
ric  and  has  identical  solutions  for  even/odd  p  (lines  I,  III  in  Fig.  2).  For  odd  n,  the 
problem  is  antimetric  and  has  identical  solutions  for  odd/even  p  (line  II  in  Fig. 2).  We 
have  choosen  k  =  2m7r  to  highlight  the  character  of  the  observation.  The  same  effect 
occurs,  however,  for  any  k  if  either  k  or  kh  are  close  to  integer  multitudes  of  tt. 

Returning  to  the  estimate  in  Theorem  3.1.  -  here:  eq  (4.1)  -  we  next  present  numer¬ 
ical  data  illuminating  the  dependence  of  the  constant  Ca  on  p.  Namely,  for  p  =  1 ...  6 
and  hk  =  1  (case  n  =  A:  =  50)  we  record  the  relative  error,  the  measured  value  of 
Ci,measip)  =  |e|i,mea*(2p)2p  and  the  predicted  C'i(p)  from  eq  (3.5). 

Table  1:  Constant  Ca{p)  in  the  Approximation  Theorem:  magnitude  as  computed 
from  measured  data  compared  to  magnitude  as  computed  from  theoretical  prediction 
in  eq  (3.5) 


p 

^  l,mea5 

Ca^meas 

Ca,  (3.5) 

1 

0.2823 

0.5646 

1.02 

2 

0.367E-1 

0.5872 

1.17 

3 

0.3095E-2 

0.6685 

1.43 

4 

0.1965E-3 

0.9049 

1.81 

5 

0.9829E-5 

0.9829 

2.33 

6 

0.4135E-6 

1.2357 

3.02 

We  see  that,  for  the  particular  case  computed,  both  the  magnitudes  of  the  measured 
constant  and  it’s  growth  rates  with  p  are  lower  than  the  upper  theoretical  estimates. 
We  conclude  that  for  the  example  under  consideration,  no  further  growth  with  p  occurs 
in  the  relative  error  due  to  the  ratio  |u|p+i/|u|i. 
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For  graphical  illustration  we  compare  in  Fig.  3  (for  k  =  SOtt  and  m  —  50)  the 
estimated  and  measured  errors  for  p  =  1 . . .  6.  Namely,  we  plot  (setting  Ca{p)  =  1) 

est„(p)  =  (n  •  2p)~^ 

for  n  =  50  and  1  <  p  <  6.  We  compare  with  the  relative  error  |e|i  as  obtained  from 
computation.  The  measured  error  is  in  close  agreement  with  the  estimator,  indicating 
that  the  estimate  (4.1)  is  sharp. 

Finally,  we  relate  the  error  of  best  approximation  to  the  number  of  degrees  of 
freedom  of  the  discrete  model.  In  one  dimension,  this  number  is  d  :=  For  fixed 

jb,  the  error  estimate  (4.1)  becomes 

Igl,  <  (4.2) 

Observe  in  Fig.  4  both  the  predicted  rate  of  convergence  w.r.  to  d  and  the  decrease  of 
the  multiplicating  factor  with  the  increase  of  p. 


Error  of  the  finite  element  solution:  The  error  of  the  finite  element  solution  is 
computed  by 


|u  ^/e|i 

E  (4'’cp'  +  aplcP  - 

j=i  1=3 


(4.3) 


where  ^ 

ej  =  |u|j  —  h  ^  (DjuDjUh  +  DjUhDjU  —  DjUhDjUh)  (4-4) 

j=i 

is  the  error  of  piecewise  linear  approximation  (cf.  part  I)  and  are  the  coefficients 
of  the  bubble  modes  in  the  local  finite  element  ansatz,  cf.  subsection  3.2. 

In  Fig.  5,  the  relative  error  of  the  finite  element  solution  is  plotted  against  the 
relative  error  of  the  best  approximation.  The  wavenumber  is  A:  =  307r  and  the  results 
are  compared  for  p  =  1 ...  6.  We  clearly  see  the  optimal  convergence  of  the  finite 
element  solution  for  sufficiently  small  h.  In  the  given  example,  the  optimality  constant 
C  in  Theorem  3.5  is  asymptotically  1  as  the  figure  shows. 

To  illustrate  the  behaviour  in  the  preasymptotic  range,  consider  the  horizontal  line 
drawn  at  e  =  0.1.  Compared  to  the  asymptotic  behavior  with  the  optimality  constant 
close  to  1,  the  finite  element  solution  is  significantly  polluted  on  the  preasymptotic 
error  level.  We  give  the  numerical  results  in  Table  2. 
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The  error  estimate  for  oscillating  solutions  has  been  given  in  Corollary  3.2,  eq  (3.98). 
For  our  solution,  the  special  case  1.  applies,  hence  we  have  the  estimate 


e  1  <  C\{p)9  +  ® 


with 


In  this  inequality,  Ci{p)6  estimates  the  error  of  best  approximation  as  discussed  in 
the  previous  paragraph.  The  second  member  C2{p)k6^  reflects  the  pollution  and  is  of 
the  same  order  as  the  phase  lag  (cf.  Theorem  3.2).  Theoretically,  the  constant  C2 
may  significantly  grow  with  p,  cf.  Theorem  3.7  vs.  Lemma  3.5.  As  the  table  shows, 
we  do  not  observe  this  growth  in  the  example  considered.  As  commented  before,  we 
cannot  exclude  that  the  theoretically  established  dependence  of  C'2(p)  on  p  is  due  to  a 
technicality  in  the  proof. 

Finally  we  observe  that  the  number  of  elements  for  which  the  finite  element  error 
is  a  fixed  magnitude  (given  tolerance)  decreases  significantly  with  increasing  p. 

Table  2:  Errors  of  finite  element  solution  and  best  approximation  at  e/e  «  0.1  for 

p  =  1 ...  6. 


®We  neglect  the  term  of  eq  (3.98) 
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p 

1 

2 

3 

4 

5 

6 

^(meas) 

^/e 

0.099947 

0.09956 

0.09683 

0.0911 

0.0829 

0.09833 

n  (S  of  elem.) 

491 

76 

35 

22 

16 

12 

of  DOF 

491 

152 

105 

88 

80 

72 

o-OfY 

0.09597 

0.09612 

0.0904 

0.0822 

0.07092 

0.7861 

^(meas) 

4a 

0.05538 

0.05607 

0.05640 

0.05543 

0.0520 

0.0626 

^fe  ^ba 

0.044567 

0.04349 

0.04040 

0.03567 

0.02380 

0.0380 

ke^ 

0.8681 

0.8707 

0.7702 

0.6373 

0.474 

0.5823 

C(p)  = 

0.0513 

0.050 

0.0525 

0.056 

0.0502 

0.0656 
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F.  Ihlenburg  and  I.  Babuska 


Figure  1:  Relative  error  (eba)  of  the  best  ^-approximation  in  S'J(n)  to  the  exact 
solution.  Rates  of  convergence  in  .H’^-seminorm  for  p  =  1, 2, . . . ,  6  (A:  =  50,  n  =  h~^  - 
number  of  elements). 


Finite  Element  Solution  to  the  Helmholtz  Equation  II ... 


Relative  error  of  best  approximation  in  Hl-seminorm,  kss  30*Pi,  n=4.70.1,  p  =  1,6,1 


Figure  2:  Relative  error  of  the  best  approximation  in  if^-seminorm  on  coarse  gr 
k  =  SOtt. 
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F.  Ihhnburg  and  1.  Babuska 


Figure  3:  Relative  error  of  the  best  approximation  in  ^^-seminorm:  estimated  vs. 
measured  values  for  k  =  SOtt,  n  =  50  and  p  =  1, 2, . . . ,  6. 
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Figure  4:  Relative  error  of  best  approximation  vs,  number  of  degrees  of  freedom  for 
fc  =  50  and  p  =  1 , 2, . . . ,  6. 
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Figure  5:  Relative  error  of  the  finite  element  solution  versus  best  approximation  error 
for  k  =  SOtt  and  p  =  1, 2, . . . ,  6. 
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